PIV_5:Dynamic model decompositio

2019-08-19 本文已影响0人闪电侠悟空

0.简要背景介绍

POD分解（PCA主成分分析）的一个显著不足是没有对时序建模的能力。热力学第二定律告诉我们时间是很重要的因素，流体力学的流动状态也是随着时间而变化的，如果PCA的基能够随着时间变化，这就太好了。
考虑时间，一阶线性动力学系统是最简单的考虑时间因素的动态系统。
DMD首先被应用到流场分析中，流场分析的天然大数据特性，也具有时间变化特性（比如湍流发展），所以是一种综合考虑时间和空间的流场分析利器（比如流场预测与控制）。由于线性动力学系统和PCA分析的可解释性，今天，大数据出现在越来越多的其他学科应用场合，DMD因此也成为分析空间复杂，时间变化的数据信号的一种新型手段。
参考资料动态模态分解(DMD)与数据科学和 DMD wiki和论文Multi-Resolution Dynamic Mode Decomposition
别人的实现代码

1248340160.jpg

1. 一阶线性动力学系统 (DMD中的Dynamic)

一阶微分方程 $\frac{d\mathbf{x}(t)}{dt}=\mathbf{K}\mathbf{x}(t)$ ,其中 $K$ 为常量。

高等数学提供了封闭解： $\mathbf{x}(t)=e^{\mathbf{K}t}\cdot \mathbf{x}(0)$ 。这个解就先放这里,后面貌似也没用，现在继续看下面一种离散的迭代解。
离散迭代表达： $\mathbf{x}(t_{i+1}) = A \cdot \mathbf{x}(t_{i})$ ，这个就是上面的微分方程的离散表达， $t_i$ 是指连续时间 $t$ 的离散采样。注意， $A$ 和 $\mathbf{K}$ 有着千丝万缕的联系，一种非正式的表达式 $\mathbf{A} = (\mathbf{I}+dt\cdot \mathbf{K})$ 。

这个 $A$ 是非常重要的，明显是一种状态转移矩阵，类比下HMM的模型 $(\lambda, A,B)$ 。直接求取和分析 $A$ 面临一个主要的问题，A太大了，比如 $100\times 100$ 的中小型流场， $\mathbf{x}(t_i)\in \mathbb{R}^{20000}$ ， $A\in \mathbb{R}^{20000\times 20000}$ 根本没法处理。直观感觉这个 $x(t_i)$ 需要降维，才能继续。

2. Snapshot的表示

为了方便起见，列向量 $\mathbf{x}(t_i)$ 可以简记为 $\mathbf{x}_i$ ，表示在 $t_i$ 时刻的拍照(snapshot)。从1到 $T$ 时刻的Snapshot当然可以用一个矩阵表示

$X_{1}^{T}=[\mathbf{x}_1,\mathbf{x}_2,\mathbf{x}_3,...,\mathbf{x}_T]$
微分方程的离散迭代表达式可以表达为矩阵形式：
$X_{2}^{T} = AX_{1}^{T-1}$
所以A可以通过广义逆矩阵（ $\cdot^ +$ ）形式化的给出
$A=X_{2}^{T}(X_{1}^{T-1})^ +$
由于前面提到的问题 $A$ 太大了，求取 $A$ 很难，对 $A$ 进行分析（比如特征值分解）也很难。

3. 对矩阵 $A$ 进行特征值分解的重要性

根据离散表达式 $\mathbf{x}_{i+1} = A \cdot \mathbf{x}_{i}$ ， $\mathbf{x}_{T} = A^T \cdot \mathbf{x}_{0}$ ，所以 $A$ 的特征值表示了对应特征向量产生的序列是收敛（ $\lambda<1$ ），发散（ $\lambda>1$ ），震荡（ $\lambda$ 虚部不为0）三种状态。特征向量就是DMD的初始时刻的投影基。所以说，得到矩阵 $A$ 的特征值分解比得到矩阵 $A$ 本身更加有用。这里也形式化的给出 $A$ 的特征值分解表达式子
$AW= W \Lambda$

4. 如何求到 $A$ 的特征值和特征向量？

$A$ 之所以难求在于 $A$ 太大，本质在于 $\mathbf{x}_i$ 的维度太高。POD分解（或者PCA，SVD分解）提供了一种降维处理的方式。先搞一个低维版本看看，说不定就帮助把高维问题给解决掉了。

$X_{1}^{T-1}=U\Sigma V^ *$ 选择这个矩阵奇异值分解，（注意只选择K阶的表达），广义逆就好处理了
$A=X_{2}^{T}(X_{1}^{T-1})^ += X_{2}^{T} V\Sigma^{-1}U^*$

以上都只是对 $X_{1}^{T-1}$ 进行了降维处理（具体的操作为 $U^*X_{1}^{T-1}$ ），但是矩阵 $A$ 还是那个 $A$ ，能不能对 $A$ 也降维处理下？所以构造出了新的降维后的 $\tilde{A}$ .
$\tilde{A} = U^*AU$

分析下矩阵大小： $\mathbf{x}_i\in \mathbb{R}^N$ ， $X_{2}^{T},X_{1}^{T-1} \in \mathbb{R}^{N\times T}$ ， $(X_{1}^{T-1})^ + \in \mathbb{R}^{T\times N}$ ， $U\in \mathbb{R}^{N\times K}$ ， $V\in \mathbb{R}^{K\times T}$ ， $\Sigma\in \mathbb{R}^{K\times K}$ ， $A\in \mathbb{R}^{N\times N}$ ， $\tilde{A}\in \mathbb{R}^{K\times K}$ ，果然实现了矩阵维度的大幅度下降，可以处理了。
分析下特征值的关系： $A$ 和 $\tilde{A}$ 的特征值是相同的。推导过程见下面。
分析下特征向量的关系：推到见下面

4.1 推导特征向量和特征值的关系

$A\mathbf{w}= \mathbf{w} \lambda$ ， $\mathbf{w}\in \mathbb{R}^N$ ，（1）
$\tilde{A}\mathbf{\tilde{w}}= \mathbf{\tilde{w}} \tilde{\lambda}$ ， $\mathbf{\tilde{w}}\in \mathbb{R}^K$ ，（2）
结论是 $\mathbf{w} =X_{2}^{T} V\Sigma^{-1} \mathbf{\tilde{w}}$ ， $\lambda=\tilde{\lambda}$ .
证明如下：令 $Z=X_{2}^{T} V\Sigma^{-1}\in \mathbb{R}^{N\times K}$ ，则 $A=ZU^*$ ， $\tilde{A}=U^*Z$ . 那么 $A$ 的特征向量也应该能够用 $Z$ 的列(行)空间张成， $\mathbf{w} = Z\mathbf{y}$ 。
结果 $AZ\mathbf{y}=ZU^*Z\mathbf{y}= Z\mathbf{y} \lambda$ , 应用 $Z$ 的左逆，得到 $\mathbf{y} = \mathbf{\tilde{w}}$ .

5. 如何用DMD分解呢？

5.1 DMD 预测与重构

回顾线性系统的离散解 $\mathbf{x}_{i+1} = A \cdot \mathbf{x}_{i}$ ，目前显然迭代式可以写成显示表达式 $\mathbf{x}_i = A^i \mathbf{x}_0$ . 同时 $\mathbf{x}_0$ 可以分别投影到基上 $\mathbf{x}_0=\Sigma_{k=0}^k<\mathbf{x}_0,\mathbf{w}_k>\mathbf{w}_k$ .
$\mathbf{x}_i = A^i \mathbf{x}_0 = A^i\Sigma_{k=0}^k<\mathbf{x}_0,\mathbf{w}_k>\mathbf{w}_k=\Sigma_{k=0}^K<\mathbf{x}_0,\mathbf{w}_k>A^i\mathbf{w}_k=\Sigma_{k=0}^K<\mathbf{x}_0,\mathbf{w}_k>\lambda_k^i\mathbf{w}_k$
简单点，并且记作矩阵形式：
$^{DMD}\mathbf{x}_i =\Sigma_{k=0}^K<\mathbf{x}_0,\mathbf{w}_k>\lambda_k^i\mathbf{w}_k=\Psi\Lambda^iB$ ，上标 $DMD$ 表示DMD重构的snapshot i。

其中 $\Psi$ 是特征向量，是处理模态的空间变化，本质和ＰＯＤ的基模态是一样的；
其中 $\Lambda^i$ 是与时间相关的量，用于处理模态的时间变化，这个是线性系统一样的；
$B$ 是和初始状态有关的量；
空间模态 $K$ 的每一个阶次都能抽离出来，单独分析。基本模态可分别展示出随着时间的发展。

5.2 python代码

时间晚了，回去吃饭，明天来写一个代码测试下

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns

sns.set()

# define time and space domains
x = np.linspace(-10, 10, 100)
t = np.linspace(0, 6 * np.pi, 80)
dt = t[2] - t[1]
Xm, Tm = np.meshgrid(x, t)

# create three spatiotemporal patterns
f1 = np.multiply(20 - 0.2 * np.power(Xm, 2), np.exp((2.3j) * Tm))
f2 = np.multiply(Xm, np.exp(0.6j * Tm))
f3 = np.multiply(5 * np.multiply(1 / np.cosh(Xm / 2), np.tanh(Xm / 2)), 2 * np.exp((0.1 + 2.8j) * Tm))

# combine signals and make data matrix
D = (f1 + f2 + f3).T

# create DMD input-output matrices
X = D[:, :-1]
Y = D[:, 1:]

# print(D)

# fig = plt.figure()
# ax = Axes3D(fig)
# ax.plot_surface(Xm, Tm, D.T, rstride=1, cstride=1, cmap='rainbow')

plt.figure()
ax = plt.axes(projection='3d')
ax.plot_surface(Xm, Tm, np.abs(D.T), cmap='rainbow')
ax.set_xlabel('x')
ax.set_ylabel('t')
plt.title("Spatial Temporal Data")
# plt.show()

# DMD核心部分 ********************************************************************
# SVD of input matrix
U2, Sig2, Vh2 = np.linalg.svd(X, False)

# rank-3 truncation
r = 3
U = U2[:, :r]
Sig = np.diag(Sig2)[:r, :r]
V = Vh2.conj().T[:, :r]

# 展示基本模态
plt.figure()
for i in range(r):
    plt.plot(U[:, i], label=str(i))
plt.legend(["1st", "2nd", "3rd"], loc='best')
plt.title("Basic POD modal")
# plt.show()

# build A tilde
Atil = np.dot(np.dot(np.dot(U.conj().T, Y), V), np.linalg.inv(Sig))
mu, W = np.linalg.eig(Atil)
plt.figure()
plt.scatter(mu.real, mu.imag)

theta = np.linspace(0, np.pi * 2, 200)
plt.plot(np.cos(theta), np.sin(theta), '--')
plt.title("Eig value")
# plt.show()

# build DMD modes (矩阵A的特征向量)
Phi = np.dot(np.dot(np.dot(Y, V), np.linalg.inv(Sig)), W)
plt.figure()
for i in range(r):
    plt.plot(Phi[:, i], label=str(i))
plt.legend(["1st", "2nd", "3rd"], loc='best')
plt.title("Eig Vector of A")
# plt.show()

# compute time evolution
b = np.dot(np.linalg.pinv(Phi), X[:, 0])  # 根据初始时刻算的B
Psi = np.zeros([r, len(t)], dtype='complex')
for i, _t in enumerate(t):
    Psi[:, i] = np.multiply(np.power(mu, _t / dt), b)

plt.figure()
plt.subplot(311)
plt.plot(Psi[0, :].real)
plt.plot(Psi[0, :].imag)
plt.subplot(312)
plt.plot(Psi[1, :].real)
plt.plot(Psi[1, :].imag)
plt.subplot(313)
plt.plot(Psi[2, :].real)
plt.plot(Psi[2, :].imag)

# plt.figure()
# for i in range(r):
#     plt.plot(Phi[:, i], label=str(i))
# plt.legend(["1st", "2nd", "3rd"], loc='best')
# plt.title("Eig Vector of A")
# # plt.show()

# compute DMD reconstruction
D2 = np.dot(Phi, Psi)
np.allclose(D, D2)  # True
plt.figure()
plt.imshow(np.abs(D2 - D))
plt.colorbar()
plt.title("Error of Reconstruction")
plt.show()

PIV_5:Dynamic model decompositio

0.简要背景介绍

1. 一阶线性动力学系统 (DMD中的Dynamic)

2. Snapshot的表示

3. 对矩阵 $A$ 进行特征值分解的重要性