Koopman-EDMD实现四旋翼非线性系统辨识与数据驱动控制

简介：本资源是一套面向控制理论与机器人方向高年级本科生及研究生的数据驱动建模与控制实践材料，聚焦四旋翼无人机这一典型非线性系统，提供基于Koopman算子与扩展动态模式分解（EDMD）的完整Matlab实现方案，解决传统建模依赖精确物理参数、非线性控制器设计复杂等实际问题。压缩包共90个文件，含41个核心Matlab函数（如edmd/eval_EDMD/get_basis/main等）、11幅结果可视化图（png/fig格式，涵盖特征值谱、MPC轨迹跟踪、EDMD误差评估等）、3个实测数据集（mat格式）及配套说明文档（md、pptx、txt），总大小48.94MB，结构清晰、模块解耦，支持2014–2024多版本Matlab运行。已有55人学习下载，代码采用参数化设计并附详尽注释，覆盖数据采集、基函数构造、Koopman矩阵辨识、模型预测控制（MPC）闭环仿真全流程，特别适合作为现代控制课程实验、毕业设计或科研入门的可复现技术模板。 Koopman-EDMD这套东西，我断断续续折腾了大半年。最早是在看非线性动力学相关文献时注意到这个词，后来发现四旋翼无人机控制圈子也在用，才意识到它不是个纯理论玩具，而是真的能把“非线性系统辨识+控制”这条路走通。这篇文章就把我在这套方法上的完整实现过程、Matlab代码思路、以及踩过的坑一次性讲清楚。

先交代一下这篇内容能干什么：如果你手里有一架四旋翼，有它的输入输出数据（或者有仿真模型可以生成数据），想通过数据驱动的方式得到一个足够准确的预测模型，然后在这个模型上做MPC、LQR之类的控制，那么Koopman-EDMD正好是适合你的路线。它核心思路很朴素——把非线性系统通过“升维观测函数”映射到一个高维空间，在这个空间里动力学近似是线性的。别被这个描述吓到，后面我会用Matlab一步步拆开给你看。

1. 从非线性到近似线性：Koopman-EDMD到底解决了什么问题

1.1 非线性控制为什么难，Koopman怎么解

四旋翼的动力学用线性系统近似只能在小角度悬停附近勉强成立。一旦你做大幅机动、快速爬升、或者带负载变化，线性化模型失配就非常明显，后面控制器设计的所有理论保证都会打折。传统处理手法无非两条路：一是精确建模，把气动阻尼、陀螺效应、电机延迟全写进方程，然后做非线性控制设计；二是反馈线性化，用微分同胚变换把非线性项消掉。这两条路都依赖“我手上这个模型是准的”，而实际中模型参数往往难测，气动系数更是随飞行状态飘。

Koopman算子的思路换了条路。它不追求在原始状态空间建模非线性，而是把状态和输入一起映射到一组观测函数张成的空间里，在这个新空间里，状态演化可以写成线性矩阵乘法的形式。做个不太严谨但很直观的类比：你在三维空间里看到的非线性轨迹，换个更高维的视角去看，可能是直线运动。Koopman-EDMD就是通过数据，自动帮你找到这个“更高维的视角”。

这点在实际工程上有很大的吸引力——如果能拿到一个近似的线性模型，那LQR、MPC、H-infinity这些非常成熟的线性控制理论就都能直接往上套，而那些理论在工业界打磨了几十年，可靠性和工程可用性都不是非线性方法能比的。

1.2 控制系统版本的EDMD推导

先定义一下我们要解决的问题。考虑离散时间的控制仿射系统：

x(k+1) = F(x(k), u(k))

其中x是系统状态，u是控制输入。传统系统辨识直接去拟合这个F，而EDMD的做法是先选定一个字典（一组标量观测函数）：

Ψ(x) = [ψ₁(x), ψ₂(x), ..., ψN(x)]ᵀ

注意这里每个ψi是一个从状态空间映射到实数的函数。常见的选法是取状态本身的各分量、二次项、交叉项、以及高次项，有时也加径向基函数。具体怎么选，我后面会专门讲，那是EDMD好坏的关键。

EDMD的核心假设是：升维之后的向量Ψ(x)在Koopman算子的作用下近似线性演化。写成公式：

Ψ(x(k+1)) ≈ K Ψ(x(k))

对于控制系统，我们要把输入也纳入。常用的处理方式是把输入增广到状态里，定义增广向量z = [x; u]，然后在每次更新时把u当成已知的外部输入，在数据中记录z(k)和z(k+1)，其中z(k)的输入部分是u(k)，z(k+1)的输入部分是u(k+1)——这里有个细节，做离线EDMD估计时，输入轨迹是已知的，所以可以直接这样构造快照对。实际部署时，会有专门的处理方式，后面第三节我会详细写。

于是我们要解决的问题变成了：找矩阵K，使得：

Ψ(z(k+1)) ≈ K Ψ(z(k))

对所有数据点成立。这就是个标准最小二乘问题。把数据堆起来——Ψ(X̃)是字典作用在增广快照上的矩阵，Ψ(Ỹ)是字典作用在下一步增广快照上的矩阵——求解：

K = Ψ(Ỹ) Ψ(X̃)ᵀ (Ψ(X̃)Ψ(X̃)ᵀ)⁻¹

在Matlab里通常写成K = G \ A的形式，其中G = Ψ(X̃)Ψ(X̃)ᵀ，A = Ψ(X̃)Ψ(Ỹ)ᵀ。这就是扩展动态模态分解（EDMD）求Koopman矩阵估计的核心。

1.3 临界澄清：Koopman线性化不等于局部线性化

这一节我必须单独提出来讲，因为很多人容易搞混。Koopman-EDMD（理论上）给出的是全局线性嵌入，也就是在观测函数空间里，动力学是精确线性的，不需要“小扰动”假设。这是它和经典局部线性化模型（比如Jacobian线性化）的本质区别。代价是观测空间的维度通常远高于原始状态维度——你以维度换取了线性度。

不过在实际数据驱动实现中，由于函数字典有限、数据有限，我们得到的K是真实Koopman算子的一个低秩近似，所以“近似线性”才是准确的说法。但即使如此，它的有效范围通常比某个工作点附近的Taylor展开大得多。我做过对比测试：同一组四旋翼模型，用Jacobian线性化得到的模型在姿态角偏差超过20度时预测就明显失真，而用Koopman模型（40维观测空间）在±50度范围内都能保持较高的预测精度。这就是升维线性化的价值所在。

2. 四旋翼无人机的动力学模型与训练数据生成

2.1 用的简化模型长什么样

很多做Koopman控制的人会直接用真实飞行数据，这没问题，但我在起步阶段强烈建议先用一个带参数的仿真模型生成数据，因为这样做有以下好处：真值可以对照、噪声水平可控、数据量可以无限扩展。等算法在仿真里跑通了，再切换到真实数据，会少很多排查问题的困扰。

我用的是经典的简化四面体模型。状态向量取12维：

x = [x, y, z, vx, vy, vz, φ, θ, ψ, p, q, r]ᵀ

控制输入取4维：

u = [T, τφ, τθ, τψ]ᵀ

其中T是总推力，后三个是滚转、俯仰、偏航力矩。动力学方程如下：

dp/dt = v dv/dt = [0, 0, -g]ᵀ + (R(φ,θ,ψ) [0,0,T]ᵀ) / m dφ/dt = p + (q sinφ + r cosφ) tanθ dθ/dt = q cosφ - r sinφ dψ/dt = (q sinφ + r cosφ) / cosθ I · dω/dt = τ - ω × (I · ω)

其中R(φ,θ,ψ)是旋转矩阵，I是惯性张量，τ = [τφ, τθ, τψ]ᵀ。注意欧拉角速率和机体角速度之间的关系是非线性的，描述的是“绕ZYX顺序旋转”下的运动学关系。这里我不展开每个公式的推导，但提醒一句：如果你的最终目标是悬停附近控制，这组方程的分量里俯仰、滚转的欧拉角速率近似等于p、q；如果是做大机动仿真，就不能这么简化。

在Matlab里我用的是ode45对连续动力学做数值积分，但训练数据和控制设计都基于离散模型。采样周期我取dt = 0.02s（50Hz），这个频率对大多数四旋翼姿态控制是合理的。如果你所用的飞行控制器是500Hz或者1kHz，建议把采样周期压缩到0.002～0.005s，但那样数据量会很大，计算EDMD矩阵时我要很小心内存。

2.2 激励信号设计：让数据“充分覆盖”

EDMD本质上是数据驱动的回归，所以训练数据的激励质量直接决定了模型好坏。关于激励，我踩过最大的坑是——只用随机噪声做激励，结果Koopman模型在小机动内预测很好，一做大机动就崩。原因很简单：随机白噪声的功率谱是平的，它很难把四旋翼的非线性动力学特征激发出来，特别是大角度姿态机动时的耦合效应。

我的做法是采用“多频正弦扫频 + 伪随机阶跃”的混合激励。具体来说：

• 姿态通道（滚转、俯仰、偏航力矩）：叠加多个不同频率的正弦信号，频率范围从0.1Hz扫到5Hz，幅值逐渐增大；
• 高度通道（总推力）：用一个带低通滤波的伪随机二元序列（PRBS），模拟起飞、降落、爬升、下降的切换过程；
• 每个通道单独扫频完成后，再叠加所有通道信号做一个综合激励，目的是让模型学习到通道间的耦合。

之所以这样做，是因为Koopman-EDMD和神经网络系统辨识类似——模型只能学会它见过的动力学区域。如果你的训练数据只覆盖了姿态角±10度的小范围，那模型在大角度下预测烂是必然的。所以在数据生成阶段，就要想清楚你的控制器今后会在什么状态空间范围内工作，然后让训练数据覆盖这个范围的1.5倍以上，给模型留出泛化余量。

%% 激励信号生成示例（伪代码） fs = 50; t = 0:1/fs:20; % 20秒数据 tau_phi = 0.02*sin(2*pi*0.5*t) + 0.015*sin(2*pi*1.3*t) + 0.01*sin(2*pi*3.7*t); tau_theta = 0.02*sin(2*pi*0.4*t + 1) + 0.012*sin(2*pi*2.1*t) + 0.008*sin(2*pi*4.2*t); T = 0.5 + 0.15*prbs_lowpass(t, 0.3); % PRBS叠加在悬停推力附近

2.3 数据预处理，直接影响矩阵条件数

EDMD的矩阵求解涉及构造G = Ψ(X)Ψ(X)ᵀ，这个矩阵的条件数直接决定了数值求解的稳定性。在实际操作中，我发现以下几项预处理几乎是必须做的，不做的话模型质量会明显下降：

第一点是去均值。如果数据不减去均值，状态向量里的常值分量（比如重力补偿后的悬停推力）会造成字典矩阵的第一列极其占优，最小二乘解会被这个偏置支配，导致动力学部分拟合不充分。我的做法是：用整个训练集的均值做中心化，然后把中心化后的数据再做EDMD。预测时，先把状态中心化，用K矩阵预测，最后再加回均值。

第二点是归一化。四旋翼状态量纲不统一：位置是米（量级可能是0.1～10），角速度是rad/s（量级可能是0.01～1），推力是N（量级可能是5～20）。如果直接把这些量混合进一个字典，数值上大的量会主导误差，模型训练时会“偏爱”学习推力通道的动力学，而忽略角速度通道。我的做法是对每个状态分量做线性归一化到[-1, 1]区间。归一化参数用训练集的均值和标准差，存下来，预测和控制阶段用同一组参数。

第三点是时间步长一致。这听起来是废话，但我真的见过有人用不同的采样率混着采集数据，最后EDMD矩阵算出来整个不对。所有数据快照对必须是同一个固定步长。

做完这三步，再去看cond(G)和rcond(G)，条件数通常能比未预处理时小几个数量级，这对后面的正则化策略选择也有很大影响。

3. EDMD在Matlab中的分步实现

3.1 字典函数库设计

这是EDMD最像“手艺活”的部分。字典选得好不好，直接决定升维后的线性模型能多精确地逼近原非线性动力学。理论上Koopman算子需要无穷维观测空间才能精确线性化，实际做EDMD只能取有限维，所以字典要尽可能捕获原系统的主要非线性特征。

我给四旋翼用的字典结构如下：

Ψ(x) = [x(1..12); % 原始状态 x(1..12).^2; % 二次项 交叉项x(i)*x(j); % 部分手动挑选的交叉项 φ^3, θ^3, ψ^3; % 姿态角三次项 vx*φ, vy*θ, vz*ψ; % 线速度和姿态的耦合项 p*φ, q*θ, r*ψ; % 角速度和姿态的耦合项 径向基函数 exp(-||x - c_k||^2/σ²)] % 几个中心点

为什么这么选？有三个理由：

• 原始状态本身一定要包含，这样才能直接读出预测的状态；
• 二次和三次项是动力学中常见的非线性（比如旋转矩阵里的sin和cos展开时会出现这些阶次）；
• 径向基函数用来补偿全局多项式拟合不到的局部非线性特征，增强升维空间的表达能力。

关于中心点c_k的选择，我用了k-means聚类，从训练数据里选8个聚类中心。这里要注意的是，径向基的中心点应该覆盖训练数据的分布范围，如果中心点设置得过于集中在某个区域，其他区域的拟合能力就会下降。

字典的总维度我一般控制在30~50维。有人可能会问，为什么不直接上100维甚至200维？我实测下来的感觉是，维度越高对训练数据的拟合越好，但过拟合风险也在增加，同时EDMD矩阵的规模变大、条件数变差，数值不稳定问题会集中爆发。我用40维左右的效果已经很理想，更多维度带来的收益边际递减。

3.2 核心代码：从快照对到Koopman矩阵

下面这是核心代码。我先构造数据矩阵X_all和Y_all，分别是当前时刻和下一时刻的状态快照集合。对于控制输入，我用增广的方式。

%% ===== EDMD核心：从观测数据估算Koopman矩阵 ===== % 输入: % X_data: nX x N 矩阵，每个列向量是当前时刻状态 % U_data: nU x N 矩阵，每个列向量是当前时刻输入 % Y_data: nX x N 矩阵，每个列向量是下一时刻状态 % dict_func: 函数句柄，返回升维观测向量 % ================================================ function [Koopman, A_koopman, B_koopman] = edmd_control(X_data, U_data, Y_data, dict_func) N = size(X_data, 2); % 构造增广快照：z = [x; u] Z = [X_data; U_data]; % 对每个时刻，计算升维观测向量 Psi_Z = zeros(dict_size, N); Psi_Y = zeros(dict_size, N); for k = 1:N Psi_Z(:, k) = dict_func(Z(:, k)); % 注意：Y_data(k+1) 对应的是下一时刻状态，但下一时刻的输入 u(k+1) 在数据生成时 % 是已知的。如果按照严格定义，Psi(Y) 应该是字典作用在 z(k+1)=[x(k+1); u(k+1)] 上。 % 许多实现里简化掉这一项，直接用 Psi([x(k+1); u(k)])，因为输入在这一步之后是待设计的。 % 我这里采用简化的处理方式，后面的控制器设计会说明为什么这样可行。 Psi_Y(:, k) = dict_func([Y_data(:, k); U_data(:, k)]); end % 计算 Gram 矩阵和交叉矩阵 G = Psi_Z * Psi_Z' / N; A = Psi_Y * Psi_Z' / N; % 求解 Koopman 矩阵（用伪逆更稳，后面解释） Koopman = A / G; % 相当于 A * pinv(G) % 分解为状态和输入两个部分（如果字典结构允许） % 这里假设字典的最后一组是输入本身，以便分离出 B 矩阵 % 具体分离方式和字典结构相关，见正文说明 A_koopman = Koopman(:, 1:nX_dict); B_koopman = Koopman(:, nX_dict+1:end); end

这里有个非常关键的理解点。标准EDMD的回归目标是：

Ψ(x(k+1)) ≈ K Ψ(x(k), u(k))

也就是说，我们用一个K矩阵把当前增广观测映射到下一步的观测。这样一旦有了K，预测时就分两步：给定当前状态和输入，构造增广观测Ψ([x;u])，乘上K得到Ψ(x_next)的估计，再从Ψ里读出状态分量。

但在控制器设计阶段，输入u是我们需要优化的变量，它不应该参与从Ψ到Ψ的静态映射，而应该被分离出来作为线性模型里的输入项。所以我的字典设计里，最后几个维度就是输入u本身（线性项）。这样K矩阵可以拆成两部分：

K = [A_koopman, B_koopman]

其中A_koopman作用在状态观测部分，B_koopman作用在输入部分。于是升维状态z = Φ(x)的演化方程写为：

z(k+1) = A_koopman * z(k) + B_koopman * u(k)

这就是可以直接丢给LQR/MPC的线性系统。这里要说明一个假设：由于我在数据回归时对Ψ([x(k+1); u(k)])做了一点简化处理（输入部分保持为u(k)而非u(k+1)），所以分离出的B矩阵隐含了“输入通过当前状态影响下一步观测”这一关系，这在大多数控制周期远大于输入变化周期的场景下是合理近似。

3.3 模型验证：单步/多步预测误差

模型估计完之后，不能急着去设计控制器，先做预测验证。这一步非常重要，也是我从反复试错中总结出来的——用EDMD模型做1步预测和做10步预测的误差表现可以差好几个量级，这直接反映了模型是否真正学到了系统动力学。

我的验证流程是这样的：

%% 模型验证 z1 = zeros(nZ, 1); z1(1:nX) = X_data(:, 1); err_total = 0; horizon = 50; % 预测50步（1秒） for startIdx = 1:100 % 取一段测试数据 x_true = X_data(:, startIdx:startIdx+horizon); u_true = U_data(:, startIdx:startIdx+horizon-1); % 初始化预测 z_pred = Phi(x_true(:,1)); x_pred = C * z_pred; % 多步预测 for k = 1:horizon z_pred = A_koopman * z_pred + B_koopman * u_true(:,k); x_pred(:, k+1) = C * z_pred; end % 计算误差 err = x_pred - x_true; err_total = err_total + sum(err(:).^2); end rmse = sqrt(err_total / (100 * horizon * nX));

这里的C矩阵是从升维观测中提取原始状态分量的投影矩阵。如果字典的前nX个维度就是原始状态本身，C就是[I_nX, 0]。

我最常遇到的情况是：单步预测RMSE很小，多步预测却发散。这通常说明模型学到了一个“瞬时正确但累积错误”的动力学——原因可能是字典维度不够、或者Koopman矩阵的谱分布与原系统不匹配。排查时可以画出Koopman矩阵的特征值分布，跟原非线性系统在典型平衡点线性化的特征值做对比。如果Koopman特征值里有离散单位圆外的点，那多步预测发散是必然的，这时要考虑加正则化或者增大数据覆盖范围。

另外要特别检查数据的边界行为。我试过一次，模型在训练数据覆盖范围内预测非常好，但稍微外推一点点（比如设定一个训练时没出现过的外部扰动），模型立刻崩掉。这是数据驱动方法逃不过的宿命——它只对你训练过的动力学区域负责。所以下游控制器设计时，要么把工作范围限制在训练覆盖区间内，要么在控制律里加保护性约束。

4. 在Koopman模型上设计数据驱动控制器

4.1 状态提升与观测重构

拿到A_koopman和B_koopman之后，你的系统变成了线性形式：

z(k+1) = A_koopman * z(k) + B_koopman * u(k)

这里的z是升维观测向量（比如40维），而实际状态x是12维。你要么直接在这个40维空间里设计控制器，然后把控制结果投影回原始状态；要么把模型降维到原始状态空间。前者是Koopman控制器的标准做法，因为线性系统的所有理论在任意维度都成立，40维对LQR/MPC来说毫无压力。

实际操作中要处理一个问题：控制器的状态反馈需要知道当前时刻的z(k)，但我们只能测到原始状态x(k)。解决办法很简单——用字典函数对测量值做提升：

z_measured = Phi(x_measured);

这相当于一个非线性观测器，不打引号的确定性变换。好在字典是显式函数，所以从x到z的变换是精确的，不需要额外设计状态观测器。如果测量噪声比较大，可以再加一个Kalman滤波器，在z空间里做状态估计，这会平滑掉一部分观测噪声对控制的影响。

4.2 离散LQR和线性MPC两种实现路线

Koopman控制器和传统线性控制器在实现上没有本质区别，核心差异只是你在用哪个矩阵控制系统。我分别试了离散LQR和线性MPC，下面把两套方案的Matlab实现要点都说一下。

方案一：离散LQR。

在升维空间里给定权重矩阵Q（针对z）和R（针对u），解离散Riccati方程：

[K_lqr, S, e] = dlqr(Ad, Bd, Q, R);

控制律就是u(k) = -K_lqr * (z_sp - z(k))，其中z_sp是参考点的升维表示。这里有个设计细节：Q矩阵不是直接在原始状态空间里设，而是要通过Q_z = C' * Q_x * C来构造（如果C = [I, 0]，意味着状态部分权重直接作用在原始状态上）。因为我们的字典前几维就是原始状态，所以这个转换很自然。

我用的权重：

Q_x = diag([10, 10, 10, 2, 2, 2, 50, 50, 20, 3, 3, 3]); % 位置、速度、姿态、角速度 R = diag([0.1, 0.1, 0.1, 0.1]); % 四个控制通道

这样位置误差的权重比角速度误差大得多，控制器会优先保证位置收敛，姿态只是中间变量。实测效果：在仿真中对阶跃参考信号的位置跟踪上升时间约0.8秒，稳态误差几乎为零，姿态角在机动过程中最大偏差约10度（在训练数据覆盖范围内）。

方案二：线性MPC。

MPC的好处是可以显式加输入约束。四旋翼的输入是有物理边界的（推力不能为负，电机转速有上下限），MPC能把这些约束直接优化进去。我用Matlab的quadprog手写了一个简化版MPC，预测时域取10步，控制时域取3步，每步计算QP问题：

function u_opt = koopman_mpc(z_current, r_ref, A, B, Q_mpc, R_mpc, u_min, u_max, Np, Nc) % 构造QP矩阵（这里用简化的分块矩阵形式） % 最小化 (Z - Z_ref)' * Q_bar * (Z - Z_ref) + U' * R_bar * U % 其中 Z = A_bar * z_current + B_bar * U [A_bar, B_bar] = build_prediction_matrices(A, B, Np, Nc); z_ref_bar = repmat(z_ref, Np, 1); H = B_bar' * Q_bar * B_bar + R_bar; f = (A_bar * z_current - z_ref_bar)' * Q_bar * B_bar; options = optimoptions('quadprog', 'Display', 'off'); U_opt = quadprog(H, f, [], [], [], [], ... repmat(u_min, Nc, 1), repmat(u_max, Nc, 1), [], options); u_opt = U_opt(1:nU); % 只取第一步 end

MPC比LQR多了约束处理能力，但代价是每步要解一个QP，实时性要求高的时候要优化求解器。我在仿真里跑50Hz的MPC，用quadprog每步大约花3~5毫秒，对非实时仿真完全够用。如果要做嵌入式实时部署，建议改用C代码生成或OSQP这类高效求解器。

4.3 闭环仿真：跟踪性能与鲁棒性对比

为了验证Koopman控制器到底行不行，我做了一组对照实验：分别用Koopman-LQR、Koopman-MPC和基于Jacobian线性化的LQR，去跟踪同一组参考轨迹，然后在系统里注入10%的参数不确定性（质量偏差和气动阻尼偏差），观察三种控制器的表现。

结果很有意思，直接列数据：

这里的轨迹是一个“绕圈+爬升+横滚机动”的组合，训练数据覆盖了这个范围。可以看到Koopman-LQR比Jacobian-LQR好不少，而MPC因为带约束优化，输入一直没饱和，所以跟踪误差最小。我觉得这个结果很有说服力——Koopman模型的优势在机动范围稍大时立刻体现出来，而Jacobian线性化在中等机动下就已经撑不住了。

鲁棒性测试更值得讲。我故意在模型中加了15%的质量误差和额外的气动阻尼力矩，Koopman-MPC还是能保持稳定，只是跟踪误差略有上升。这说明Koopman模型学到的动力学结构具有一定的泛化能力，不完全依赖精确参数。这个特性和神经网络控制器有点像，但Koopman模型的输出是显式线性矩阵，天然继承了线性控制理论的可解释性和稳定性分析方法，这点是神经网络给不了的。

5. 实测中的坑和调参经验

5.1 基函数选的不是越多越好

我一开始做EDMD时，抱着“多加点基函数总是好的”的想法，把字典加到80多维，什么七次项、各种交叉项统统塞进去。结果G矩阵条件数直接到1e15以上，模型预测一塌糊涂，数值上全是噪声。

后来才明白，字典函数的选取要匹配系统的本质非线性结构，而不是盲目堆高阶项。四旋翼动力学里的非线性主要来自旋转矩阵（三角函数）和刚体动力学中的叉乘项（角速度耦合），这些非线性在物理上可以用低次多项式合理近似。我最终采用的字典是高次项只保留到三次，交叉项也只加有物理意义的组合（比如角速度×姿态角）。加完跑出来rcond(G)从1e-16提升到了1e-6左右，模型预测立刻正常了。

我的经验是：在字典维度增加之前，先去看你正在处理的非线性结构长什么样，再决定加什么基函数。如果不知道怎么选，可以先从多项式阶次递增试起，每次加一个阶次看验证误差的变化，找到误差开始反弹的点，那就是过拟合的起点。

5.2 持续激励与数据量，怎么平衡

EDMD本质上做的是最小二乘回归，所以数据必须满足持续激励条件——也就是说，升维后的观测矩阵Ψ(X)的行要是满秩的，这样G矩阵才是可逆的。如果激励不充分，你会发现cond(G)很大、K矩阵的元素数值很奇怪、预测发散。

数据量方面，我用下来发现不是越多越好，但太少了一定不行。比较合理的经验是：数据长度N至少要是字典维度N_Ψ的5~10倍。比如40维字典，至少需要200~400个快照对（按0.02s采样就是4~8秒数据）。实际中我总是生成3~5倍于最低要求的数据量，然后做随机打乱，避免时间相关性影响回归效果。

还有一个细节很容易忽略：数据集里的状态分布是否均匀。如果你让无人机大部分时间悬停在一个固定点，只有偶尔几次机动，那训练数据里悬停附近的样本会占绝大多数，EDMD的误差会被那些“平凡样本”稀释，真正有信息量的机动样本反而被淹没。我的做法是去掉训练集中过于相似的样本（用简单的空间网格化采样），或者对训练样本做加权，让不同状态区域的样本贡献平衡。

5.3 数值稳定性：正则化、伪逆与条件数

EDMD计算K矩阵时，我见过很多人的代码直接写K = A * inv(G)，这在条件数差一点的情况下就是个灾难。inv是显式求逆，数值上极不稳定。正确做法是用伪逆或者mldivide：

K = A / G; % 相当于 A * pinv(G)，但数值更稳 % 或者 K = A * pinv(G, tol); % 显式指定容差

进一步的，我在实际实现中加入Tikhonov正则化：在求解时给G加一个小的对角扰动：

K = A * pinv(G + lambda * eye(size(G)));

这个lambda要调，一般取1e-6到1e-3之间。这样做能显著降低K矩阵的方差，尤其是当数据量不够或者字典维度偏高时。代价是引入微小偏差，但实测下来偏差在可接受范围内，换来的是预测稳定性大幅提升。

另外建议用rcond函数诊断矩阵健康状况。我在代码里做了这样一个检查：

if rcond(G) < 1e-10 warning('Gram矩阵接近奇异，请检查数据激励或增加正则化'); end

这个检查能帮你在早期发现数据问题，而不是等到模型预测一团糟才回头排查。

5.4 我的一些实际体会和后续扩展思路

Koopman-EDMD这条技术路线，我一年用下来的总体感觉是：它不适合当作“银弹”去解决所有的非线性控制问题，但在“中等非线性强度+数据丰富”的场景里，它的性价比非常高。相比神经网络黑箱控制器，它的模型形式是可解释的、可分析的、可以嫁接所有线性控制理论；相比精确机理建模，它又不需要你花几个月去辨识参数。

如果你想在这个基础上继续延展，我建议关注几个方向：

一是做在线自适应。把EDMD改成递归形式（递归EDMD），每收到一批新数据就更新K矩阵，让模型跟随系统退化或环境变化。这个在无人机长期飞行场景下很有价值。

二是结合鲁棒控制。Koopman模型有模型误差，可以在设计控制器时把估计误差的界算进去，然后做鲁棒MPC（Tube-MPC），这样即使模型不完美，闭环也有理论保证。

三是用实验数据替代仿真数据。如果飞行场地和硬件条件允许，用真实飞行日志训练Koopman模型。这时需要格外注意数据的时间对齐、传感器噪声建模和激励信号的可实现性，这些都会直接影响模型的最终效果。

我自己接下来打算把这套框架接到一个嵌套的级联结构上——内环用传统PID做姿态控制，外环用Koopman-MPC做位置控制。原因是姿态通道的动力学非线性更强、响应更快，用传统方法已经非常成熟；而位置通道的强耦合和大范围机动更适合Koopman模型来提升精度。这种分层混合的思路，在很多实际工程问题里往往是性价比最高的方案。

本文还有配套的精品资源，点击获取