DMD (Dynamic Mode Decomposition)

DMD分解

DMD（动态模态分解）是一种数据驱动的算法，可以用来处理时间序列，用这些时间序列数据构造一个近似的线性动态系统模型，在一定条件下逼近或预测非线性系统的行为，并可以预测未来的行为。

Standard DMD

数据预处理和 POD 算法基本一致（但不必减去平均流场） POD (Proper Orthogonal Decomposition) 即有一组时间快照，时间间隔为 $\varDelta t$

$\begin{matrix} x_1& x_2& \cdots& x_{n+1}\\ \end{matrix}$

并假设这组数据可以由线性动态系统产生，即

$x_{k+1}=Ax_k$

这个方程类似差分方程求解，具体可以看求解差分方程 | 哈基启的博客

现在目标是要找到矩阵 $A_{m\times m}$ 并求解特征向量和特征值

将这组数据分成两组矩阵

$X_{m\times n}=\left[ \begin{matrix} x_1& x_2& \cdots& x_n\\ \end{matrix} \right] \,\,\,\,\,\,\,\,\,\,\,\,X_{m\times n}^{'}=\left[ \begin{matrix} x_2& x_3& \cdots& x_{n+1}\\ \end{matrix} \right]$

对 $X$ 进行 reduced SVD 分解，r 为矩阵 $X$ 的秩

$X_{m\times n}=U_{m\times r}\varSigma _{r\times r}V^T_{r\times n}$

所以矩阵 $A$ 就可以如下表示

$A=X'X^+=X'V\varSigma ^{-1}U^T$

其中 $X^+=V\varSigma ^{-1}U^T$ 时矩阵 $X$ 的伪逆

定义 $A$ 的相似矩阵 $\widetilde{A}_{r\times r}$

$\widetilde{A}=U^TAU=U^TX'V\varSigma ^{-1}$

进行 reduced SVD 分解，并构造矩阵 $A_{m\times m}$ 的相似矩阵 $\widetilde{A}_{r\times r}$ ，是为了在后续求解特征值和特征向量时更加方便，因为原矩阵 $A_{m\times m}$ 的维度通常较大

计算矩阵 $\widetilde{A}_{r\times r}$ 的特征值和特征向量

$\widetilde{A}w=\lambda w$

并对 $\widetilde{A}_{r\times r}$ 进行 EVD 分解

$\widetilde{A}W_{r\times r}=W\varLambda_{r\times r}$

矩阵 $A$ 的特征向量（也是 DMD 的模态 mode）为

$\hat{\phi }=Uw$

此时定义的 DMD 的模态叫做 projected DMD modes

Exact DMD

现考虑一组数据对 $\left\{ \left( x_1,y_1 \right) ,\left( x_2,y_2 \right) ,\cdot \cdot \cdot ,\left( x_n,y_n \right) \right\}$ ，将其分成两组矩阵

$X_{m\times n}=\left[ \begin{matrix} x_1& x_2& \cdots& x_n\\ \end{matrix} \right] \,\,\,\,\,\,\,\,\,\,\,\,Y_{m\times n}=\left[ \begin{matrix} y_1& y_2& \cdots& y_n\\ \end{matrix} \right]$

相比于 Standard DMD，Exact DMD 在原始数据上弱化了时序性，只强调 $x_i,y_i$ 是成对出现

并且 Standard DMD 可以看成 Exact DMD 的特殊情况，即 $y_k=z_k\ \ \ x_k=z_{k-1}$

$Y=AX$

现在目标是要找到矩阵 $A_{m\times m}$ 并求解特征向量和特征值

剩下的步骤和 1.1 节完全相同，唯一不同的是最后一步定义 DMD 模态

$\phi =\frac{1}{\lambda}YV\varSigma ^{-1}w$

也有的定义为

$\phi =YV\varSigma ^{-1}w$

两个都表示同一个模态，只是前面的系数不同，此时的模态被称为 Exact DMD modes

证明： $\phi =YV\varSigma ^{-1}w$ 是 $A$ 的特征向量（前面系数不会影响特征向量所以只用证一个）

$X^+=V\varSigma ^{-1}U^T$

$A=YX^+=YV\varSigma ^{-1}U^T=BU^T$

其中 $B=YV\varSigma ^{-1}$

$\widetilde{A}=U^TYV\varSigma ^{-1}=U^TB$

设 $\widetilde{A}w =\lambda w$ ，让 $\phi=Bw$

$A\phi =BU^TBw=B\widetilde{A}w=\lambda Bw=\lambda \phi$

两种模态的对比

projected DMD modes

$\hat{\phi }=Uw$

Exact DMD modes （现在一般都采用这种定义）

$\phi =YV\varSigma ^{-1}w$

$\mathbb{P}_X$ 是投影到 $X$ 的列空间上的投影矩阵， $\hat{\phi }$ 是矩阵 $\mathbb{P}_XA$ 的特征向量，对应的特征值是 $\lambda$ ，并且 $\hat{\phi }=\mathbb{P}_X\phi$

$X_{m\times n}=U_{m\times r}\varSigma _{r\times r}V^T_{r\times n}$ 矩阵 $X$ 进行了 SVD 分解， $U^TU$ 是 $r\times r$ 的单位矩阵，矩阵 $\mathbb{P}_X=UU^T$ 是投影矩阵，投影到 $U$ 的列空间上，又因为 $U_{m\times r}$ 是 C(X) 的一组标准正交基，所以矩阵 $\mathbb{P}_X=UU^T$ 也是投影到 $X$ 的列空间上的投影矩阵。

可以通过最小二乘法理解投影矩阵：

设方程组 $Ux=b$ 无解，可以通过最小二乘法找到 $b$ 在 $U$ 列空间上的投影

$Ux=b$

$U^TU\hat{x}=U^Tb$

$I\hat{x}=U^Tb$

$\hat{x}=U^Tb$

向量 $b$ 在 $U$ 的列空间上的投影就是向量 $P$

$P=U\hat{x}=UU^Tb$

可以看出 $UU^T$ 是投影矩阵

$\mathbb{P}_XA\hat{\phi}=\left( UU^T \right) \left( BU^T \right) \left( Uw \right) =U\left( U^TB \right) w=U\widetilde{A}w=\lambda Uw=\lambda \hat{\phi}$

所以 $\hat{\phi }$ 是矩阵 $\mathbb{P}_XA$ 的特征向量

并且

$U^T\phi =U^TBw=\widetilde{A}w=\lambda w$

$\mathbb{P}_X\phi =UU^T\phi =U\lambda w=\lambda \hat{\phi}$

由上式可以看出模态 $\hat{\phi}$ 是模态 $\phi$ 在 $X$ 列空间上的投影，且有 $U^T\phi=\lambda w$ , $U^T\hat{\phi}=w$

又由矩阵 $X$ 的 POD 分解， $U_{m\times r}$ 也是 $X$ 的前 r 个 POD 模态，所以也可以理解为模态 $\hat{\phi}$ 是模态 $\phi$ 在前 r 个 POD 模态上的投影。

但注意 DMD 的模态没有正交性

结果整理

现在得到了 DMD 模态 $\phi$ 和对应的特征值 $\lambda$ ，对应的模态的矩阵为 $\varPhi_{m\times r}$ ，特征值写成对角阵 $\varLambda_{r\times r}$

原动态系统为

$x_{k+1}=Ax_k$

若 $r=m$ 即 $\varPhi_{m\times m}$ 且 $\varPhi$ 可逆，则该动态系统的解可写为

$x_k=\varPhi \varLambda ^{k-1}\varPhi ^{-1}x_1$

但是我们只找到了前 r 个特征向量（模态），所以

$x_k=\varPhi \varLambda ^{k-1}\mathbf{b}$

我们需要寻找向量 $b_{r\times 1}$ ，这个向量可以理解为初始条件在 DMD 模态上的最优投影。

它可以通过以下最小二乘问题求得：

$\mathbf{b} = \arg\min_{b} \|x_1 - \Phi \mathbf{b}\|_2$

最终

$x_k=\varPhi \varLambda ^{k-1}\mathbf{b}=\underset{j=1}{\overset{r}{\varSigma}}\phi _j\lambda _{j}^{k-1}b_j$

其中 $\phi _j$ 是 DMD 模态（矩阵 $A$ 的特征向量）， $\lambda_{j}$ 是 DMD 特征值（矩阵 $A$ 的特征值）， $b_j$ 是模态幅值（mode amplitude）

上面的 DMD 模型是以离散时间步长 $\varDelta t$ 为基础的，即描述动态系统在快照之间的离散演化规律：

$x_{k+1}=Ax_k$

若用连续时间微分方程形式来表示系统的演化：

$\frac{dx}{dt}=Bx$

其解为 $x(t)=e^{Bt}x_1=\varPhi e^{\varOmega t}\varPhi^{-1}x_1=\varPhi e^{\varOmega t}\mathbf{b}$ 求解微分方程 | 哈基启的博客

其在特征值上的关系为 $\omega =\frac{\ln \left( \lambda \right)}{\varDelta t}$

我们就能将离散时间的 DMD 模型转化为连续时间的动力系统形式，即从

$x_k=\varPhi \varLambda ^{k-1}\mathbf{b}=\underset{j=1}{\overset{r}{\varSigma}}\phi _j\lambda _{j}^{k-1}b_j$

转换为

$x\left( t \right) =\varPhi e^{\varOmega t}\mathbf{b}=\underset{j=1}{\overset{r}{\varSigma}}\phi _je^{\omega _jt}b_j$

定义：

第 $i$ 个模态的增长率 $g_i$ 和频率 $f_i$ 为

$g_i=\frac{\text{Re}\left\{ \ln \left( \lambda _i \right) \right\}}{\varDelta t}\ \ \ \ \ \ \ f_i=\frac{\text{Im}\left\{ \ln \left( \lambda _i \right) \right\}}{2\pi \varDelta t}$

其中 $\text{Im}\left\{ \ln \left( \lambda _i \right) \right\} =arg\left( \lambda _i \right)$

最终时间快照序列

$X=\left[ \begin{matrix} x_1& x_2& \cdots& x_n\\ \end{matrix} \right]$

可以写为

$X=\varPhi D_bV_{and}=\left[ \begin{matrix} \phi _1& \phi _2& \cdots& \phi _r\\ \end{matrix} \right] \left[ \begin{matrix} b_1& & & \\ & b_2& & \\ & & \ddots& \\ & & & b_r\\ \end{matrix} \right] \left[ \begin{matrix} 1& \lambda _1& \cdots& \lambda _{1}^{n-1}\\ \vdots& \vdots& & \vdots\\ \vdots& \vdots& & \vdots\\ 1& \lambda _r& \cdots& \lambda _{r}^{n-1}\\ \end{matrix} \right]$