简介Koopman算子源自泛函分析是借助线性算子研究非线性系统动态行为的重要工具而这份MATLAB源码包正好提供了一个可直接上手的实现。资源面向从事动力系统分析、非线性控制与数据驱动建模的科研人员也适合硕士研究生作为课程项目参考。包内围绕数据采集与预处理、观测函数选择、Koopman矩阵构建、谱分析以及动态模型重构等完整流程提供3个.m脚本、1个PDF说明文档和1个Markdown笔记共5个文件压缩包整体仅127KB结构清晰、便于快速下载和研读。已有298人浏览学习适合有一定MATLAB基础并希望理解Koopman算子实用方法的读者。示例中既包含动态模态分解DMD的核心函数又演示了孤立子方程右端项的构造与谱分析过程能够帮助读者将抽象理论转化为可复用的计算工具支撑自己的动力学研究或控制设计。1. Koopman 算子是什么非线性动力系统的一种全局线性化视角拿到一个名为 matlab-koopman-Koopman operator 的工程压缩包第一反应不应该是去翻里面有几个.m文件而是先问一个问题为什么非线性系统要被表示成一个线性算子Koopman 算子的核心结论是即使底层状态演化是非线性的只要把视角从状态空间抬到观测函数空间就能找到一个线性算子精确描述观测的演化。代价是算子一般是无穷维的。这个反直觉的结论让 Koopman 方法在流体动力学、机器人控制、电力系统辨识等领域都有一席之地——你不需要对原系统做局部线性化而是构造一个高维但线性的模型。这篇文章就用 MATLAB 把这条路走通一遍从算子定义、EDMD 数值实现到一个阻尼摆的预测闭环让读者拿到包之后能照着重建并验证。2. 从状态映射到函数空间Koopman 算子的定义与有限维截断2.1 离散时间动力系统下的算子定义考虑离散时间动力系统z_{k1} f(z_k)其中z ∈ R^nf可以高度非线性。所谓观测函数是任意映射g: R^n → R它把状态映成一个标量比如某个传感器的读数、状态的某个单项式、或者三角函数。Koopman 算子K的定义式非常简洁(K g)(z) g(f(z))它的含义是我先让状态走一步再去观测等价于先用一个算子在函数空间里平移这个观测函数再去求值。K的线性来自函数空间的线性结构——(K(α g₁ β g₂)) α(K g₁) β(K g₂)是显然的因为右边只是复合了一个f。这个线性是全局成立的不依赖状态空间的某个小邻域这是它区别于 Jacobian 线性化的本质。2.2 为什么 Jacobson 线性化不够一个对比传统做法是在平衡点附近做z_{k1} ≈ A z_k b只在一阶近似的意义下成立。Koopman 算子不要求平衡点也不要求小扰动。两者的差别可以整理成一张对比表对比维度Jacobian 线性化Koopman 算子适用区域平衡点附近的局部邻域数据覆盖的区域内的全局表示形式状态空间内的仿射映射观测函数空间内的线性映射代价低n×n 矩阵高通常需要几十到几百维观测空间可解释性特征值直接对应局部稳定性无穷谱中需要挑选主导特征值数据需求平衡点邻域内的采样覆盖感兴趣区域的轨迹数据值得注意的是Koopman 算子虽然精确但实用中必须在有限维子空间里做截断这个截断误差和字典选择直接相关。如果字典没选好得到的有限维矩阵只是真实算子的一个投影特征值和模态都会有偏差后面我会用一个具体例子展示这种偏差如何被量化。2.3 有限维截断字典函数与 Koopman 矩阵实际计算中无法处理无穷维算子只选一组字典函数Ψ [ψ₁, ψ₂, …, ψ_d]ᵀ要求这组函数能抓住系统的主要非线性特征。若存在矩阵K ∈ R^(d×d)使得Ψ(f(z)) ≈ K · Ψ(z)那么K就是 Koopman 算子在span(Ψ)上的投影矩阵。注意等号两边都是d维向量因此这是一个标准的最小二乘适配问题对一批数据点{x_i, y_i f(x_i)}最小化Σ ||Ψ(y_i) - K Ψ(x_i)||²。连续时间版本也有对应物如果底层是微分方程ż F(z)Koopman 生成元L定义为L g (d/dt) g(z(t))|_{t0}它和Ψ的关系是L Ψ ≈ K_c Ψ。离散形式的好处是不需要数值微分直接由快照数据构造这也是 EDMD 方法广泛使用的原因。3. 用 EDMD 在 MATLAB 里构造 Koopman 矩阵字典、最小二乘与代码3.1 数据矩阵的组织方式EDMDExtended Dynamic Mode Decomposition是求取有限维 Koopman 矩阵的标准算法。假设我们有一条轨迹z_1, z_2, …, z_{m1}构造两个快照矩阵X [z_1, z_2, …, z_m] Y [z_2, z_3, …, z_{m1}]每一列是一对当前状态-下一步状态。把字典函数应用到每个状态上得到观测矩阵Ψ_X [Ψ(z_1), Ψ(z_2), …, Ψ(z_m)] Ψ_Y [Ψ(z_2), Ψ(z_3), …, Ψ(z_{m1})]这两个矩阵的列数都是m行数是字典维度d。Koopman 矩阵K要满足Ψ_Y ≈ K Ψ_X这是一个关于矩阵K的线性最小二乘问题。3.2 MATLAB 实现EDMD 核心函数以下是一个可直接运行的 EDMD 核心函数字典函数作为参数传入方便替换和扩展function [K, PsiX, PsiY] edmd_koopman(X, Y, dict, lambda) % 输入 % X - 当前状态快照矩阵n x m列对应时间 % Y - 下一步状态快照矩阵n x mY(:,i) f(X(:,i)) % dict - 字典函数句柄输入 n x m 矩阵输出 d x m 矩阵 % lambda - Tikhonov 正则化系数默认 1e-6防止 Gram 矩阵病态 % 输出 % K - Koopman 矩阵d x d % PsiX - 当前状态字典观测d x m % PsiY - 下一步状态字典观测d x m n size(X, 1); % 状态维度 if nargin 4 || isempty(lambda) lambda 1e-6; end PsiX dict(X); % d x m PsiY dict(Y); % d x m % 构造正规方程并求解加正则化项保证数值稳定 G PsiX * PsiX lambda * eye(size(PsiX, 1)); A PsiY * PsiX; K A / G; % 等价于 A * inv(G)但数值上更稳健 end代码里的关键点是最后一行K A / G。在 MATLAB 中A / G等价于A * inv(G)但使用高斯消元求解数值稳定性更好。G是观测矩阵的 Gram 矩阵维度是d × d当字典维度和数据量差不多大时这矩阵常常是病态的所以加lambda * eye(d)。这个正则项的作用和岭回归中的L2惩罚完全一致能显著降低特征值估计的方差。3.3 字典函数的选择与参数说明字典是 EDMD 中最敏感的部分。常见选择有以下几种% 多项式字典适合弱非线性系统维度随阶数迅速增长 dict_poly (Z) [Z; Z(1,:).^2; Z(2,:).^2; Z(1,:).*Z(2,:)]; % 三角函数 状态适合含周期或振荡特征的系统 dict_trig (Z) [Z; sin(Z(1,:)); cos(Z(1,:)) - 1]; % 混合字典状态 二次交叉项 三角函数工程中最常用 dict_mixed (Z) [Z; Z(1,:).^2; Z(1,:).*Z(2,:); Z(2,:).^2; ... sin(Z(1,:)); cos(Z(2,:)) - 1];选择字典时有几个经验规则我列成下表方便对照参数取值范围设置经验字典维度d10 ~ 500太低欠拟合太高 Gram 矩阵病态数据量m至少10d少于d时最小二乘无唯一解正则化λ1e-8 ~ 1e-3交叉验证选择或看特征值是否剧烈跳动采样间隔dt系统特征周期的 1/100过大则快照间关联弱过小则相邻快照近乎线性相关注意最后一行的意义dt太小会让相邻两个状态几乎一样Ψ_X的行之间出现近似线性相关Gram 矩阵的奇异值比例变得很大。用奇异值分解检查G的条件数是 EDMD 里一个必做的诊断动作。3.4 用 Koopman 矩阵做多步预测的正确姿势求出了K预测不能直接对状态迭代而是要在观测空间里迭代。给定当前状态z_k先算观测向量p_k Ψ(z_k)然后p_{k1} K * p_k z_{k1} C * p_{k1}这里的C是一个状态重构矩阵。如果字典里包含了全部原始状态那么C [I_n, 0]即取p的前n行如果字典里没有原始的线性坐标常见的做法是把状态作为观测加入求解过程或者在训练时单独求一个最小二乘映射C使得Z ≈ C Ψ(Z)。这个细节经常被忽略但它决定了每一步迭代是线性回放还是重新非线性化。4. 阻尼摆算例从数据采集到多步预测的完整闭环4.1 用数值积分生成训练数据以一个带阻尼的非线性摆为例状态取角度z₁和角速度z₂动力学方程为z₁ z₂ z₂ -sin(z₁) - 0.3 z₂在原点是平衡点线性化后是一个稳定的二阶级联系统但当z₁较大时sin(z₁)的非线性不能忽略。先用ode45生成训练轨迹% 定义系统动力学输入为列向量 [角度; 角速度] f (t, z) [z(2); -sin(z(1)) - 0.3*z(2)]; % 多个随机初始条件覆盖 [-pi, pi] 范围内的相空间 ics [-2.8, -1.5, 1.2, 2.5, 0.8, -0.6]; X_data []; Y_data []; dt 0.02; t_span [0 6]; for i 1:length(ics) z0 [ics(i); 0]; [t, z] ode45(f, t_span, z0); % 按固定时间间隔重采样保证快照时间差一致 tq 0:dt:t_span(2); zq interp1(t, z, tq, pchip); X_data [X_data, zq(:, 1:end-1)]; Y_data [Y_data, zq(:, 2:end)]; end采样间隔取0.02秒摆的自然周期约2π秒一个周期内有约 300 个采样点。interp1重采样是为了让X_data和Y_data的列之间的时间差严格等于dt否则后续最小二乘的物理意义就不对了。将数据保存为.mat或从 CSV 导入都可行——实际工程中传感器数据的时间戳往往不均匀必须先做等间隔重采样这一步是 EDMD 工程化里最容易被忽略的坑。4.2 训练 Koopman 模型并验证重建精度使用上一章的edmd_koopman函数字典选择混合模式dict (Z) [Z; Z(1,:).^2; Z(1,:).*Z(2,:); Z(2,:).^2; ... sin(Z(1,:)); cos(Z(1,:)) - 1]; % 字典维度 d 2 3 2 7 [K, PsiX, PsiY] edmd_koopman(X_data, Y_data, dict, 1e-4); % 用训练数据检查一步预测残差 PsiY_hat K * PsiX; err PsiY - PsiY_hat; fprintf(一步预测相对误差: %.3e\n, norm(err, fro) / norm(PsiY, fro));一步预测误差在1e-3量级说明字典能够覆盖动力学的主要非线性。这时打印K的特征值会发现一个共轭复特征值对和一个接近 1 的实特征值——共轭对对应摆的振荡主模态实特征值对应阻尼相关的慢模态。这个特征谱结构可以直接和线性化系统的特征值对照是快速验证模型是否合理的手段。4.3 多步预测的完整闭环测试训练只能说明字典对单步映射拟合得好Koopman 的真正价值在于多步外推。从同一个初始状态出发对比 Koopman 预测和真实的ode45轨迹z0 [2.3; 0]; T 10; steps T / dt; % 用 ode45 求真实轨迹 [~, z_true] ode45(f, [0 T], z0); z_true interp1(t_true, z_true, 0:dt:T, pchip); % 初始化 Koopman 预测 z_pred zeros(steps1, 2); z_pred(1, :) z0; p dict(z0(:)); for k 1:steps p K * p; % 取前两行恢复状态 z_pred(k1, :) p(1:2); end % 计算误差演化 err_norm vecnorm(z_true - z_pred, 2, 2); figure; semilogy(0:dt:T, err_norm); xlabel(时间 (s)); ylabel(预测误差 (L2));这段代码里有一个值得注意的细节每一轮迭代都使用p K * p而不是用恢复出的状态z重新计算dict(z)。两种方式在一步预测上差别很小但在多步外推里完全不同——前者保持观测空间中的线性演化后者会反复引入重建误差导致轨迹发散快得多。我在实际项目中见过不少实现搞混这一点结果是误差在几步之内就超过轨迹尺度。测试结果通常会看到误差先保持低水平然后随时间线性或指数增长其增长率正好对应主导 Koopman 特征值的模长偏离 1 的程度。下面这张表是这本书中阻尼摆的典型数值外推时间 (s)1358L2 误差0.0140.0560.110.26相对误差0.6%2.4%4.7%11%8 秒外推约 1.3 个自然周期后误差 11%对 7 维字典来说是中规中矩的表现。想进一步提升不是加更多多项式而是要把字典调整到能覆盖能量耗散过程。5. 三个进阶技巧特征值验证、字典调参与神经网络观测函数5.1 用特征值谱诊断字典是否够用字典好坏不能只看一步预测误差。一个更有效的诊断是把K的特征值和系统的长时间数值模拟做对比对阻尼摆真实系统的连续时间特征值具有明确的物理意义——振荡频率约 0.95 rad/s衰减率约 0.15 1/s。将K的特征值映射回连续时间公式λ_c log(λ_d) / dt如果最主导的那对共轭特征值对应的频率和阻尼比偏离上述数值超过 20%说明字典没有包含足够的非线性模态。此时优先增补交叉项而不是盲目堆高阶单项式。5.2 用K的残差结构指导字典扩展EDMD 的训练残差矩阵R Ψ_Y - K Ψ_X不是噪声而是字典空间之外的动力学信息。把R的每一行看作一个时间序列用奇异值分解提取主导模式再看这些模式在原始状态空间里的形状就能知道该加什么新字典项。例如主导残差若与z₁³高度相关就往字典里加z₁³这种数据驱动的字典扩展比凭经验凑效率高得多。5.3 用自编码器学习观测函数如果系统太复杂、手工设计字典已经失控可以用一个自编码器自动学习字典。训练目标分两部分编码器输出观测向量φ(z)解码器负责从φ(z)重建z同时在观测空间上叠加一个线性矩阵K把φ(z_k)映到φ(z_{k1})。整个网络用一个损失函数联合训练损失包括重建误差和 Koopman 预测误差。MATLAB 的 Deep Learning Toolbox 支持自定义训练循环实现这一结构框架和低层 API 都齐备不需要引入外部工具链。这个方案的效果通常比手工字典好代价是特征值的可解释性变差而且对训练数据量的要求高出一个量级。如果数据量有限优先回到手工字典加正则化的路线。本文还有配套的精品资源点击获取