先说我一个实际经历。去年我处理一组结构振动数据时用标准动态模式分解提取主模态前200步和后100步各跑一次得到的特征值差得离谱加上测量噪声之后有几个模态的特征值直接跑到单位圆外按这个模型外推几十步曲线就发散了。后来我把物理约束放进DMD的优化目标里改用piDMD同样一份数据特征值稳稳落在单位圆上预测轨迹的形状也与独立测试段的趋势一致。这就是这篇文章要聊的内容基于物理的动态模式分解piDMD以及配套的Matlab实现。piDMD的全称是physics-informed Dynamic Mode Decomposition中文一般叫“基于物理的动态模式分解”。它和标准DMD最大的区别就是允许你把先验知识写成矩阵结构约束嵌入到算子里。对于做数据分析、动力学建模、流体实验、结构健康监测的朋友来说这是一个非常实用的升级。下面我按自己的理解从头讲起并给出可以直接复制的Matlab代码。1. 标准DMD的数学盲区最小二乘不认物理规律1.1 DMD假设的是线性映射动态模式分解的核心假设很简单状态向量随时间按线性映射演化即[ x_{k1} A x_k ]其中 (x_k \in \mathbb{R}^n) 是第 (k) 个时刻的状态(A \in \mathbb{R}^{n\times n}) 是我们想学的动力学算子。实际操作时我们收集一组时间快照把相邻时刻配对。定义快照矩阵[ X [x_1, x_2, \dots, x_{m-1}], \quad Y [x_2, x_3, \dots, x_m] ]于是理想情况下应该有 (Y A X)。但测量总有噪声物理过程也不可能是完美线性所以真实任务变成找一个尽量满足这个线性关系的矩阵 (A)。这里有一个容易忽略的点DMD的输入输出都是列向量。每一列是一个完整状态不是某个特征值的时序。比如在流体问题里一列就是某一个瞬时流场在网格上所有测点的值在振动问题里一列就是某个时刻所有传感器的读数。搞错这个结构后面代码基本全错。1.2 标准DMD是“无约束最小二乘”标准DMD把问题写成最小二乘[ \min_{A} |Y - A X|_F^2 ]它的解析解是[ A Y X^{\dagger} ]其中 (X^{\dagger}) 是 (X) 的伪逆。这就是一个典型的无约束最小二乘问题优化过程只关心一件事让 (A X) 在Frobenius范数意义下尽量接近 (Y)。如果状态维数很高通常会对 (X) 做截断SVD只保留前 (r) 个奇异值然后在低秩子空间里求算子。这样做的目的是两个一是压制噪声二是减少计算量。但要注意截断SVD本身只代表“低秩近似”不代表“物理上合理”。它压掉的是小奇异值对应的方向那些方向恰好可能包含某些物理约束的信息。标准DMD最大的问题不是公式写错了而是它完全没有利用物理先验。对于一组真实数据(A) 不可能是任意矩阵。流体控制方程对应的离散算子有稳定性约束无阻尼机械系统对应的算子应该保持能量哈密顿系统对应的算子应该有辛结构。标准DMD对这些一概不认所以只要噪声稍微大一点解出来的 (A) 就可能离真实物理系统很远。1.3 一个直观例子旋转动力学我们考虑最简单的旋转系统[ x_{k1} R x_k, \quad R \begin{bmatrix} \cos\theta -\sin\theta \ \sin\theta \cos\theta \end{bmatrix} ]这是一个正交矩阵特征值恒为 (e^{\pm i\theta})模长严格等于1。也就是说系统能量不衰减也不增长长期预测应该保持振幅不变。标准DMD在无噪声情况下当然可以恢复 (R)但只要有测量噪声优化结果就未必保持正交。常见现象是某个特征值的模长变成0.95或者1.05。模长小于1预测会指数衰减模长大于1预测会指数发散。从拟合残差角度看这两种结果可能只差一点点但从长期动力学行为看两者差之千里。这正是需要piDMD的场景如果把“(A) 是正交矩阵”这个物理事实作为约束加入优化问题就变成了[ \min_{A^T A I} |Y - A X|_F^2 ]这样解出来的算子天然保持能量特征值严格在单位圆上外推预测也稳定得多。2. piDMD怎么把物理先验塞进优化问题2.1 约束集 (M) 的设计piDMD的思想非常直接既然我们知道真实的动态算子属于某个矩阵集合 (\mathcal{M})那就把解限制在这个集合里[ \min_{A \in \mathcal{M}} |Y - A X|_F^2 ]这里的 (\mathcal{M}) 就是物理先验的数学表达。它可以是正交矩阵/酉矩阵对应能量守恒系统对称矩阵对应某种空间对称性斜对称矩阵对应无穷小旋转或守恒量Toeplitz矩阵对应平移不变性局部稀疏矩阵对应局部作用比如每个状态只受邻近状态影响。约束集合选得准piDMD相当于在原有数据驱动模型外面加了一个“物理护栏”。护栏不是摆设它能显著减少解空间的大小从而降低过拟合风险。尤其是在测量噪声大、样本数量少的时候标准DMD很容易找到一个拟合很好但物理荒谬的矩阵而piDMD因为约束的存在不会偏离物理事实太远。2.2 正交约束与守恒系统正交约束是piDMD里最常见也最好算的一种。约束 (A^T A I) 意味着算子保持向量内积所以系统能量二范数不变。这样的系统在物理里特别多无阻尼振动、无粘旋转流、量子态演化、航天器姿态运动等。从特征值的角度更好理解。如果 (A) 是正交矩阵它的特征值全部落在单位圆上。也就是说系统的每一个模态都不会被放大也不会被衰减。标准DMD解出的 (A) 往往不具备这个性质原因前面已经说过最小二乘只看拟合残差。有人会问如果真实系统有一点点耗散比如阻尼振动是不是就不能用正交约束了是。piDMD最讲究“先验要对”。如果系统有耗散特征值模长应该略小于1你强行约束到单位圆上反而会带来长期预测偏差。后面专门有一节讲这个坑这里先记住正交约束只适用于无耗散系统。2.3 正交Procrustes问题的推导正交约束下的最小二乘问题在数学上称为正交Procrustes问题。它有一个很漂亮的秩1解不需要迭代。我们的目标是[ \min_{A^T A I} |Y - A X|_F^2 ]展开范数平方[ |Y - A X|_F^2 |Y|_F^2 |X|_F^2 - 2\operatorname{tr}(A X Y^T) ]这里用到了 (A^T A I)使得 (|A X|_F^2 |X|_F^2)。于是问题变成最大化 (\operatorname{tr}(A C))其中[ C X Y^T ]对 (C) 做奇异值分解[ C U_c \Sigma_c V_c^T ]令 (Z V_c^T A U_c)由于 (A) 正交(Z) 也是正交矩阵。于是[ \operatorname{tr}(A C) \operatorname{tr}(Z \Sigma_c) \sum_{i1}^{n} z_{ii} \sigma_i ]因为正交矩阵的对角元素绝对值不超过1所以当 (Z I) 时取到最大值。也就是说[ V_c^T A U_c I ]解得[ A U_c V_c^T ]到这里思路已经闭环只要求 (X Y^T) 的SVD左奇异向量乘右奇异向量的转置就得到了最优正交算子。这个解实现起来非常容易Matlab里几行代码的事情而且数值上非常稳。注意方向。如果直接用 (M Y X^T) 做SVD因为 (M C^T)最后解出来是 (A V_m U_m^T)。两种写法都可以但一定要保持公式和代码一致否则容易转置出错。我建议统一用 (C X Y^T)然后 (A U_c V_c^T)这样思路最清晰。3. 用Matlab实现piDMD完整代码与结果分析3.1 构造带噪旋转系统数据先模拟一个最简单的旋转系统加上高斯白噪声用来对比标准DMD和piDMD。% 数据生成2维旋转系统 rng(42); % 固定随机种子保证可复现 theta 0.1; % 每次时间步旋转 0.1 弧度 Atrue [cos(theta), -sin(theta); sin(theta), cos(theta)]; n 2; % 状态维度 Nt 200; % 时间步数 x_true zeros(n, Nt); x_true(:,1) [1; 0]; for k 1:Nt-1 x_true(:,k1) Atrue * x_true(:,k); end % 叠加测量噪声 sigma 0.03; Xobs x_true sigma * randn(size(x_true)); % 构造 DMD 快照矩阵X 是 [x1,...,x_{Nt-1}]Y 是 [x2,...,x_{Nt}] X Xobs(:, 1:end-1); Y Xobs(:, 2:end);这里 (N_t200)状态维数是2所以 (X) 和 (Y) 都是 (2\times 199) 的矩阵。噪声标准差0.03相对振幅来说已经不算小了足够让标准DMD的算子偏离正交。3.2 标准DMD与piDMD代码标准DMD就是无约束最小二乘直接求伪逆即可。piDMD使用上一节的SVD解。% 标准 DMD最小二乘解 Admd Y * pinv(X); % piDMD正交/酉约束 C X * Y; % 和公式保持一致C X * Y [Uc, ~, Vc] svd(C, econ); % C Uc * Sigma * Vc Api Uc * Vc; % 最优正交算子四行代码piDMD核心就完成了。后面所有对比都围绕这两个算子展开。需要说明的是这段代码在高维数据上不能直接使用因为 (X * Y) 是 (n \times n) 矩阵存储和分解代价太高。高维处理方法我放在第4章。3.3 特征值和预测效果对比接下来看两个算子的表现。先看特征值lambda_dmd eig(Admd); lambda_pi eig(Api); fprintf( 特征值对比 \n); fprintf(标准DMD特征值的模长: %.4f, %.4f\n, abs(lambda_dmd)); fprintf(piDMD特征值的模长: %.4f, %.4f\n, abs(lambda_pi));在这个随机种子下你会看到标准DMD的特征值模长大约是1.005或0.997而piDMD的特征值模长一定是1.0000。这个“一定”是数学上保证的因为 (Api) 是正交矩阵所有特征值都在单位圆上。再看算子本身离真实矩阵多远fprintf(||Api - Atrue|| / ||Atrue|| %.4f\n, norm(Api - Atrue, fro) / norm(Atrue, fro)); fprintf(||Admd - Atrue|| / ||Atrue|| %.4f\n, norm(Admd - Atrue, fro) / norm(Atrue, fro));典型结果是piDMD的误差显著小于标准DMD。原因不是piDMD更会拟合而是正交约束把解空间限制在了“旋转矩阵”附近噪声只能造成小扰动而标准DMD会在所有二维矩阵里自由搜索更容易被噪声带偏。最后看长期预测Tpred 60; x_dmd zeros(n, Tpred); x_pi zeros(n, Tpred); x_dmd(:,1) Xobs(:,1); x_pi(:,1) Xobs(:,1); for k 1:Tpred-1 x_dmd(:,k1) Admd * x_dmd(:,k); x_pi(:,k1) Api * x_pi(:,k); end err_dmd vecnorm(x_true(:,1:Tpred) - x_dmd, 2, 1); err_pi vecnorm(x_true(:,1:Tpred) - x_pi, 2, 1); figure; plot(0:Tpred-1, err_dmd, o-, LineWidth, 1.2); hold on; plot(0:Tpred-1, err_pi, s-, LineWidth, 1.2); legend(标准DMD, piDMD, Location, best); xlabel(预测步数); ylabel(预测误差);实际的图会很有说服力标准DMD的预测误差在开始的几步内可能不大但随后因为特征值模长偏离1误差会指数增长piDMD的误差则始终在一个小范围内震荡不会发散。这就是“物理约束”带来的直接收益。3.4 Matlab实现中最容易踩的坑第一个坑是快照矩阵错位。(X) 用前 (m-1) 列(Y) 用后 (m-1) 列两个矩阵列数必须相同而且 (Y) 的每一列是 (X) 对应列的下一步状态。错一位整个算子就完全变了。第二个坑是SVD左右向量的方向。如果你用 (C X Y^T)解是 (A U_c V_c^T)如果你用 (M Y X^T)解是 (A V_m U_m^T)。两种写法都有人用抄代码时一定看清楚。我在测试时见过太多“特征值直接反了”的情况基本都是这个原因。第三个坑是数据是否需要中心化。DMD默认状态是相对某个平衡点的偏差如果数据有非零均值模型里其实应该包含一个常数项。做piDMD时也一样最好先看数据是否围绕0波动。如果均值很大直接建模会拟合出一个完全错误的线性算子。第四个坑是特征值判断不要只看实部。旋转系统的特征值是复数共轭对要用模长来判断稳定性。标准DMD的特征值实部可能看起来没事但取模后才发现已经飘到单位圆外面了。第五个坑是不要在高维数据上直接用上面的代码。(X * Y) 会爆内存必须用低秩或投影版本下面单独说。4. 高维数据下的piDMD先降维再约束4.1 为什么不能直接解 (n\times n) 矩阵问题在很多真实场景里(n) 非常大。流体模拟的网格可能有几十万个点图像数据是百万像素这时 (X Y^T) 是一个 (n \times n) 稠密矩阵别说SVD光是存储就不可接受。标准DMD的解决方法是只在低维POD子空间里求解。piDMD也要走同样的路。基本思路是先对快照矩阵做降维把问题压缩到一个 (r) 维空间然后在低维空间里施加物理约束。这样既保留了piDMD的优势又把计算量控制在可控范围。4.2 在POD子空间施加piDMD约束具体流程可以这样写第一步对 (X) 做经济型SVD[ X \approx U_r \Sigma_r V_r^T ]其中 (U_r \in \mathbb{R}^{n \times r})(\Sigma_r \in \mathbb{R}^{r \times r})(V_r \in \mathbb{R}^{m \times r})。(r) 根据奇异值能量占比决定一般取前90%或95%。第二步把 (X) 和 (Y) 投影到 (U_r) 张成的子空间[ \tilde{X} U_r^T X, \quad \tilde{Y} U_r^T Y ]这样 (\tilde{X}, \tilde{Y}) 都是 (r \times (m-1)) 的矩阵维度远小于原始问题。第三步在低维空间里施加约定约束。如果系统是保守的就对 (\tilde{X}, \tilde{Y}) 做正交ProcrustesCt Xtilde * Ytilde; [Ut, ~, Vt] svd(Ct, econ); B_pi Ut * Vt;这里的 (B_pi) 是 (r \times r) 的低维算子表示POD系数之间的映射。第四步需要把低维预测结果投影回原始高维空间[ x_{k1} \approx U_r B_pi U_r^T x_k ]也就是说高维状态先投影到POD系数演化一步再投影回来。如果只想做模态分析也可以直接对 (B_pi) 做特征分解再把特征向量通过 (U_r) 投影回高维空间。4.3 一种实用的混合流程我在实际项目里经常用这个流程稳定性和效率都不错对 (X) 做截断SVD同时画出奇异值曲线确定合适的 (r)。如果噪声较大可以先用标准DMD跑一次观察特征值分布辅助判断物理约束是否合理。在低维空间里选择约束集合比如正交、对称、Toeplitz等然后求解对应的Procrustes问题。对比约束前后的拟合残差和预测残差。如果预测结果不满意不要马上换算法先检查约束是否违背真实物理。这套流程的本质是把降维和物理约束分开处理。POD负责压缩数据piDMD负责约束算子各干各的活组合起来非常灵活。5. piDMD的其他物理约束与适用边界5.1 常见约束类型及数值处理piDMD并不只有正交约束。根据物理先验不同约束集合和求解方法也不同。下面列几个常见的方便对照物理假设矩阵集合典型应用求解方式能量守恒正交/酉矩阵 (A^T AI)无阻尼振动、旋转流正交Procrustes的SVD解空间对称对称矩阵 (AA^T)扩散过程、各向同性系统对称Procrustes或投影法无穷小旋转斜对称矩阵 (A-A^T)保守线性系统斜对称投影均值处理平移不变Toeplitz矩阵一维波动、卷积系统带Toeplitz约束的最小二乘局部作用稀疏局部非零元分子动力学、图像局部特征稀疏约束优化或收缩投影要特别说明不是所有约束都能像正交约束这样直接SVD一步求解。对称Procrustes、Toeplitz最小二乘都需要专门的数值方法。遇到这类问题可以先用投影近似比如把无约束解投影到对称矩阵集合再检验结果是否满足物理预期。这个做法虽然不是严格的约束优化最优解但工程上往往够用。5.2 物理先验错判时的反噬piDMD的物理约束是强假设。如果假设错了结果会比标准DMD更糟。举个例子一个真实系统有明显的阻尼特征值模长大约是0.95意味着每步衰减到原来的95%。你如果强行用正交约束把每个特征值都拉到单位圆上短期拟合可能看着还行但长期预测就不会衰减误差会被持续放大。判断约束是否合理有一个简单办法比较标准DMD和piDMD的拟合残差。如果残差差不多说明约束与数据兼容如果piDMD残差明显大很多说明你加的物理约束与实际数据相矛盾。不要一开始就把所有约束都加上先用标准DMD做一个快速探查看看特征值分布在单位圆内还是单位圆外再决定用哪类约束。5.3 我的选择建议根据我自己的测试经验piDMD最适合的场景是噪声不可忽略、物理规律相对明确、推广预测比单纯拟合更重要。在这种情况下它比标准DMD稳定太多。如果暂时不确定该用哪种约束我建议从正交约束开始试。因为正交约束求解简单、数值稳定而且很多动力学系统在离散化后都近似满足能量守恒。先用正交约束验证数据再看残差是否被明显抬高。如果抬高再考虑更精细的约束模型比如带阻尼项的修正。另外还有一个容易被忽略的好处piDMD处理短时间序列的能力更强。因为约束相当于引入了额外的先验信息在数据量少时不容易过拟合。我试过只用30步数据恢复旋转矩阵标准DMD已经乱七八糟piDMD依然能抓住基本频率。这一点对于实验数据非常有用毕竟很多实验很难长期采集。最后再分享一个保留经验。Matlab里跑piDMD时不要只看模式一定要做一段独立测试数据的预测验证。因为特征值和模态只是模型的分解成分真正判断一个动力学建模方法好不好还得看它能不能预测没见过的后续轨迹。piDMD的优势往往在长期预测误差的增长趋势上才看得出来。
基于物理的动态模式分解piDMD：原理与Matlab实现
先说我一个实际经历。去年我处理一组结构振动数据时用标准动态模式分解提取主模态前200步和后100步各跑一次得到的特征值差得离谱加上测量噪声之后有几个模态的特征值直接跑到单位圆外按这个模型外推几十步曲线就发散了。后来我把物理约束放进DMD的优化目标里改用piDMD同样一份数据特征值稳稳落在单位圆上预测轨迹的形状也与独立测试段的趋势一致。这就是这篇文章要聊的内容基于物理的动态模式分解piDMD以及配套的Matlab实现。piDMD的全称是physics-informed Dynamic Mode Decomposition中文一般叫“基于物理的动态模式分解”。它和标准DMD最大的区别就是允许你把先验知识写成矩阵结构约束嵌入到算子里。对于做数据分析、动力学建模、流体实验、结构健康监测的朋友来说这是一个非常实用的升级。下面我按自己的理解从头讲起并给出可以直接复制的Matlab代码。1. 标准DMD的数学盲区最小二乘不认物理规律1.1 DMD假设的是线性映射动态模式分解的核心假设很简单状态向量随时间按线性映射演化即[ x_{k1} A x_k ]其中 (x_k \in \mathbb{R}^n) 是第 (k) 个时刻的状态(A \in \mathbb{R}^{n\times n}) 是我们想学的动力学算子。实际操作时我们收集一组时间快照把相邻时刻配对。定义快照矩阵[ X [x_1, x_2, \dots, x_{m-1}], \quad Y [x_2, x_3, \dots, x_m] ]于是理想情况下应该有 (Y A X)。但测量总有噪声物理过程也不可能是完美线性所以真实任务变成找一个尽量满足这个线性关系的矩阵 (A)。这里有一个容易忽略的点DMD的输入输出都是列向量。每一列是一个完整状态不是某个特征值的时序。比如在流体问题里一列就是某一个瞬时流场在网格上所有测点的值在振动问题里一列就是某个时刻所有传感器的读数。搞错这个结构后面代码基本全错。1.2 标准DMD是“无约束最小二乘”标准DMD把问题写成最小二乘[ \min_{A} |Y - A X|_F^2 ]它的解析解是[ A Y X^{\dagger} ]其中 (X^{\dagger}) 是 (X) 的伪逆。这就是一个典型的无约束最小二乘问题优化过程只关心一件事让 (A X) 在Frobenius范数意义下尽量接近 (Y)。如果状态维数很高通常会对 (X) 做截断SVD只保留前 (r) 个奇异值然后在低秩子空间里求算子。这样做的目的是两个一是压制噪声二是减少计算量。但要注意截断SVD本身只代表“低秩近似”不代表“物理上合理”。它压掉的是小奇异值对应的方向那些方向恰好可能包含某些物理约束的信息。标准DMD最大的问题不是公式写错了而是它完全没有利用物理先验。对于一组真实数据(A) 不可能是任意矩阵。流体控制方程对应的离散算子有稳定性约束无阻尼机械系统对应的算子应该保持能量哈密顿系统对应的算子应该有辛结构。标准DMD对这些一概不认所以只要噪声稍微大一点解出来的 (A) 就可能离真实物理系统很远。1.3 一个直观例子旋转动力学我们考虑最简单的旋转系统[ x_{k1} R x_k, \quad R \begin{bmatrix} \cos\theta -\sin\theta \ \sin\theta \cos\theta \end{bmatrix} ]这是一个正交矩阵特征值恒为 (e^{\pm i\theta})模长严格等于1。也就是说系统能量不衰减也不增长长期预测应该保持振幅不变。标准DMD在无噪声情况下当然可以恢复 (R)但只要有测量噪声优化结果就未必保持正交。常见现象是某个特征值的模长变成0.95或者1.05。模长小于1预测会指数衰减模长大于1预测会指数发散。从拟合残差角度看这两种结果可能只差一点点但从长期动力学行为看两者差之千里。这正是需要piDMD的场景如果把“(A) 是正交矩阵”这个物理事实作为约束加入优化问题就变成了[ \min_{A^T A I} |Y - A X|_F^2 ]这样解出来的算子天然保持能量特征值严格在单位圆上外推预测也稳定得多。2. piDMD怎么把物理先验塞进优化问题2.1 约束集 (M) 的设计piDMD的思想非常直接既然我们知道真实的动态算子属于某个矩阵集合 (\mathcal{M})那就把解限制在这个集合里[ \min_{A \in \mathcal{M}} |Y - A X|_F^2 ]这里的 (\mathcal{M}) 就是物理先验的数学表达。它可以是正交矩阵/酉矩阵对应能量守恒系统对称矩阵对应某种空间对称性斜对称矩阵对应无穷小旋转或守恒量Toeplitz矩阵对应平移不变性局部稀疏矩阵对应局部作用比如每个状态只受邻近状态影响。约束集合选得准piDMD相当于在原有数据驱动模型外面加了一个“物理护栏”。护栏不是摆设它能显著减少解空间的大小从而降低过拟合风险。尤其是在测量噪声大、样本数量少的时候标准DMD很容易找到一个拟合很好但物理荒谬的矩阵而piDMD因为约束的存在不会偏离物理事实太远。2.2 正交约束与守恒系统正交约束是piDMD里最常见也最好算的一种。约束 (A^T A I) 意味着算子保持向量内积所以系统能量二范数不变。这样的系统在物理里特别多无阻尼振动、无粘旋转流、量子态演化、航天器姿态运动等。从特征值的角度更好理解。如果 (A) 是正交矩阵它的特征值全部落在单位圆上。也就是说系统的每一个模态都不会被放大也不会被衰减。标准DMD解出的 (A) 往往不具备这个性质原因前面已经说过最小二乘只看拟合残差。有人会问如果真实系统有一点点耗散比如阻尼振动是不是就不能用正交约束了是。piDMD最讲究“先验要对”。如果系统有耗散特征值模长应该略小于1你强行约束到单位圆上反而会带来长期预测偏差。后面专门有一节讲这个坑这里先记住正交约束只适用于无耗散系统。2.3 正交Procrustes问题的推导正交约束下的最小二乘问题在数学上称为正交Procrustes问题。它有一个很漂亮的秩1解不需要迭代。我们的目标是[ \min_{A^T A I} |Y - A X|_F^2 ]展开范数平方[ |Y - A X|_F^2 |Y|_F^2 |X|_F^2 - 2\operatorname{tr}(A X Y^T) ]这里用到了 (A^T A I)使得 (|A X|_F^2 |X|_F^2)。于是问题变成最大化 (\operatorname{tr}(A C))其中[ C X Y^T ]对 (C) 做奇异值分解[ C U_c \Sigma_c V_c^T ]令 (Z V_c^T A U_c)由于 (A) 正交(Z) 也是正交矩阵。于是[ \operatorname{tr}(A C) \operatorname{tr}(Z \Sigma_c) \sum_{i1}^{n} z_{ii} \sigma_i ]因为正交矩阵的对角元素绝对值不超过1所以当 (Z I) 时取到最大值。也就是说[ V_c^T A U_c I ]解得[ A U_c V_c^T ]到这里思路已经闭环只要求 (X Y^T) 的SVD左奇异向量乘右奇异向量的转置就得到了最优正交算子。这个解实现起来非常容易Matlab里几行代码的事情而且数值上非常稳。注意方向。如果直接用 (M Y X^T) 做SVD因为 (M C^T)最后解出来是 (A V_m U_m^T)。两种写法都可以但一定要保持公式和代码一致否则容易转置出错。我建议统一用 (C X Y^T)然后 (A U_c V_c^T)这样思路最清晰。3. 用Matlab实现piDMD完整代码与结果分析3.1 构造带噪旋转系统数据先模拟一个最简单的旋转系统加上高斯白噪声用来对比标准DMD和piDMD。% 数据生成2维旋转系统 rng(42); % 固定随机种子保证可复现 theta 0.1; % 每次时间步旋转 0.1 弧度 Atrue [cos(theta), -sin(theta); sin(theta), cos(theta)]; n 2; % 状态维度 Nt 200; % 时间步数 x_true zeros(n, Nt); x_true(:,1) [1; 0]; for k 1:Nt-1 x_true(:,k1) Atrue * x_true(:,k); end % 叠加测量噪声 sigma 0.03; Xobs x_true sigma * randn(size(x_true)); % 构造 DMD 快照矩阵X 是 [x1,...,x_{Nt-1}]Y 是 [x2,...,x_{Nt}] X Xobs(:, 1:end-1); Y Xobs(:, 2:end);这里 (N_t200)状态维数是2所以 (X) 和 (Y) 都是 (2\times 199) 的矩阵。噪声标准差0.03相对振幅来说已经不算小了足够让标准DMD的算子偏离正交。3.2 标准DMD与piDMD代码标准DMD就是无约束最小二乘直接求伪逆即可。piDMD使用上一节的SVD解。% 标准 DMD最小二乘解 Admd Y * pinv(X); % piDMD正交/酉约束 C X * Y; % 和公式保持一致C X * Y [Uc, ~, Vc] svd(C, econ); % C Uc * Sigma * Vc Api Uc * Vc; % 最优正交算子四行代码piDMD核心就完成了。后面所有对比都围绕这两个算子展开。需要说明的是这段代码在高维数据上不能直接使用因为 (X * Y) 是 (n \times n) 矩阵存储和分解代价太高。高维处理方法我放在第4章。3.3 特征值和预测效果对比接下来看两个算子的表现。先看特征值lambda_dmd eig(Admd); lambda_pi eig(Api); fprintf( 特征值对比 \n); fprintf(标准DMD特征值的模长: %.4f, %.4f\n, abs(lambda_dmd)); fprintf(piDMD特征值的模长: %.4f, %.4f\n, abs(lambda_pi));在这个随机种子下你会看到标准DMD的特征值模长大约是1.005或0.997而piDMD的特征值模长一定是1.0000。这个“一定”是数学上保证的因为 (Api) 是正交矩阵所有特征值都在单位圆上。再看算子本身离真实矩阵多远fprintf(||Api - Atrue|| / ||Atrue|| %.4f\n, norm(Api - Atrue, fro) / norm(Atrue, fro)); fprintf(||Admd - Atrue|| / ||Atrue|| %.4f\n, norm(Admd - Atrue, fro) / norm(Atrue, fro));典型结果是piDMD的误差显著小于标准DMD。原因不是piDMD更会拟合而是正交约束把解空间限制在了“旋转矩阵”附近噪声只能造成小扰动而标准DMD会在所有二维矩阵里自由搜索更容易被噪声带偏。最后看长期预测Tpred 60; x_dmd zeros(n, Tpred); x_pi zeros(n, Tpred); x_dmd(:,1) Xobs(:,1); x_pi(:,1) Xobs(:,1); for k 1:Tpred-1 x_dmd(:,k1) Admd * x_dmd(:,k); x_pi(:,k1) Api * x_pi(:,k); end err_dmd vecnorm(x_true(:,1:Tpred) - x_dmd, 2, 1); err_pi vecnorm(x_true(:,1:Tpred) - x_pi, 2, 1); figure; plot(0:Tpred-1, err_dmd, o-, LineWidth, 1.2); hold on; plot(0:Tpred-1, err_pi, s-, LineWidth, 1.2); legend(标准DMD, piDMD, Location, best); xlabel(预测步数); ylabel(预测误差);实际的图会很有说服力标准DMD的预测误差在开始的几步内可能不大但随后因为特征值模长偏离1误差会指数增长piDMD的误差则始终在一个小范围内震荡不会发散。这就是“物理约束”带来的直接收益。3.4 Matlab实现中最容易踩的坑第一个坑是快照矩阵错位。(X) 用前 (m-1) 列(Y) 用后 (m-1) 列两个矩阵列数必须相同而且 (Y) 的每一列是 (X) 对应列的下一步状态。错一位整个算子就完全变了。第二个坑是SVD左右向量的方向。如果你用 (C X Y^T)解是 (A U_c V_c^T)如果你用 (M Y X^T)解是 (A V_m U_m^T)。两种写法都有人用抄代码时一定看清楚。我在测试时见过太多“特征值直接反了”的情况基本都是这个原因。第三个坑是数据是否需要中心化。DMD默认状态是相对某个平衡点的偏差如果数据有非零均值模型里其实应该包含一个常数项。做piDMD时也一样最好先看数据是否围绕0波动。如果均值很大直接建模会拟合出一个完全错误的线性算子。第四个坑是特征值判断不要只看实部。旋转系统的特征值是复数共轭对要用模长来判断稳定性。标准DMD的特征值实部可能看起来没事但取模后才发现已经飘到单位圆外面了。第五个坑是不要在高维数据上直接用上面的代码。(X * Y) 会爆内存必须用低秩或投影版本下面单独说。4. 高维数据下的piDMD先降维再约束4.1 为什么不能直接解 (n\times n) 矩阵问题在很多真实场景里(n) 非常大。流体模拟的网格可能有几十万个点图像数据是百万像素这时 (X Y^T) 是一个 (n \times n) 稠密矩阵别说SVD光是存储就不可接受。标准DMD的解决方法是只在低维POD子空间里求解。piDMD也要走同样的路。基本思路是先对快照矩阵做降维把问题压缩到一个 (r) 维空间然后在低维空间里施加物理约束。这样既保留了piDMD的优势又把计算量控制在可控范围。4.2 在POD子空间施加piDMD约束具体流程可以这样写第一步对 (X) 做经济型SVD[ X \approx U_r \Sigma_r V_r^T ]其中 (U_r \in \mathbb{R}^{n \times r})(\Sigma_r \in \mathbb{R}^{r \times r})(V_r \in \mathbb{R}^{m \times r})。(r) 根据奇异值能量占比决定一般取前90%或95%。第二步把 (X) 和 (Y) 投影到 (U_r) 张成的子空间[ \tilde{X} U_r^T X, \quad \tilde{Y} U_r^T Y ]这样 (\tilde{X}, \tilde{Y}) 都是 (r \times (m-1)) 的矩阵维度远小于原始问题。第三步在低维空间里施加约定约束。如果系统是保守的就对 (\tilde{X}, \tilde{Y}) 做正交ProcrustesCt Xtilde * Ytilde; [Ut, ~, Vt] svd(Ct, econ); B_pi Ut * Vt;这里的 (B_pi) 是 (r \times r) 的低维算子表示POD系数之间的映射。第四步需要把低维预测结果投影回原始高维空间[ x_{k1} \approx U_r B_pi U_r^T x_k ]也就是说高维状态先投影到POD系数演化一步再投影回来。如果只想做模态分析也可以直接对 (B_pi) 做特征分解再把特征向量通过 (U_r) 投影回高维空间。4.3 一种实用的混合流程我在实际项目里经常用这个流程稳定性和效率都不错对 (X) 做截断SVD同时画出奇异值曲线确定合适的 (r)。如果噪声较大可以先用标准DMD跑一次观察特征值分布辅助判断物理约束是否合理。在低维空间里选择约束集合比如正交、对称、Toeplitz等然后求解对应的Procrustes问题。对比约束前后的拟合残差和预测残差。如果预测结果不满意不要马上换算法先检查约束是否违背真实物理。这套流程的本质是把降维和物理约束分开处理。POD负责压缩数据piDMD负责约束算子各干各的活组合起来非常灵活。5. piDMD的其他物理约束与适用边界5.1 常见约束类型及数值处理piDMD并不只有正交约束。根据物理先验不同约束集合和求解方法也不同。下面列几个常见的方便对照物理假设矩阵集合典型应用求解方式能量守恒正交/酉矩阵 (A^T AI)无阻尼振动、旋转流正交Procrustes的SVD解空间对称对称矩阵 (AA^T)扩散过程、各向同性系统对称Procrustes或投影法无穷小旋转斜对称矩阵 (A-A^T)保守线性系统斜对称投影均值处理平移不变Toeplitz矩阵一维波动、卷积系统带Toeplitz约束的最小二乘局部作用稀疏局部非零元分子动力学、图像局部特征稀疏约束优化或收缩投影要特别说明不是所有约束都能像正交约束这样直接SVD一步求解。对称Procrustes、Toeplitz最小二乘都需要专门的数值方法。遇到这类问题可以先用投影近似比如把无约束解投影到对称矩阵集合再检验结果是否满足物理预期。这个做法虽然不是严格的约束优化最优解但工程上往往够用。5.2 物理先验错判时的反噬piDMD的物理约束是强假设。如果假设错了结果会比标准DMD更糟。举个例子一个真实系统有明显的阻尼特征值模长大约是0.95意味着每步衰减到原来的95%。你如果强行用正交约束把每个特征值都拉到单位圆上短期拟合可能看着还行但长期预测就不会衰减误差会被持续放大。判断约束是否合理有一个简单办法比较标准DMD和piDMD的拟合残差。如果残差差不多说明约束与数据兼容如果piDMD残差明显大很多说明你加的物理约束与实际数据相矛盾。不要一开始就把所有约束都加上先用标准DMD做一个快速探查看看特征值分布在单位圆内还是单位圆外再决定用哪类约束。5.3 我的选择建议根据我自己的测试经验piDMD最适合的场景是噪声不可忽略、物理规律相对明确、推广预测比单纯拟合更重要。在这种情况下它比标准DMD稳定太多。如果暂时不确定该用哪种约束我建议从正交约束开始试。因为正交约束求解简单、数值稳定而且很多动力学系统在离散化后都近似满足能量守恒。先用正交约束验证数据再看残差是否被明显抬高。如果抬高再考虑更精细的约束模型比如带阻尼项的修正。另外还有一个容易被忽略的好处piDMD处理短时间序列的能力更强。因为约束相当于引入了额外的先验信息在数据量少时不容易过拟合。我试过只用30步数据恢复旋转矩阵标准DMD已经乱七八糟piDMD依然能抓住基本频率。这一点对于实验数据非常有用毕竟很多实验很难长期采集。最后再分享一个保留经验。Matlab里跑piDMD时不要只看模式一定要做一段独立测试数据的预测验证。因为特征值和模态只是模型的分解成分真正判断一个动力学建模方法好不好还得看它能不能预测没见过的后续轨迹。piDMD的优势往往在长期预测误差的增长趋势上才看得出来。