企业官网建设流程全解析

简介：基于Koopman算子与扩展动态模式分解（EDMD）的四旋翼无人机数据驱动控制Matlab实现方案，面向计算机、电子工程及数学专业学生，可用于课程项目、综合实验与毕业课题。压缩包共包含90个文件，其中核心为41个脚本与函数文件，另配有可视化图片、仿真数据、幻灯片说明及备份内容，整体约49.66MB，并兼容Matlab 2014、2019a与2024a等版本，便于多环境使用。当前已有110人学习下载。整个方案重点演示了基于Koopman算子将非线性系统转化为线性描述的方法，以及利用EDMD从实际测量数据中提取动态模式的具体流程；代码中包含模型预测控制、EDMD评估与轨迹生成等模块，均采用参数化设计并附有详尽注释，实验数据可直接加载复现，便于理解关键数据结构与算法逻辑。通过学习，读者能够掌握四旋翼数据驱动建模与控制的核心步骤，理解特征值分析、控制性能对比等难点，为后续研究或毕业设计提供实践基础。

1. 直接给四旋翼建模为什么越建越虚：Koopman算子换了个玩法

做四旋翼无人机控制的人大概都有同感：动力学看着简单，真要把模型建到能上机的精度，质量、转动惯量、气动系数每一个都够折腾。小角度假设能把姿态方程拉成线性，但一旦做大速度机动或重载飞行，线性模型又立刻失效。这两年我在Matlab里试了一套基于Koopman算子与EDMD算法的数据驱动控制方案：绕开逐项参数辨识，直接用飞行状态数据把非线性动力学提升成高维线性模型，再在这个线性模型上设计MPC。这套方案尤其适合做控制算法预研、仿真验证和课题起步，不需要额外硬件，不依赖精确气动参数，一台装好Matlab的电脑就能把整条链路跑通。

2. 先把Koopman算子讲透：非线性系统在什么条件下能变成线性模型

2.1 Koopman算子到底在做什么：状态提升与线性演化

传统建模思路是直接写状态方程：x_{k+1} = F(x_k)。四旋翼的F包含推力与姿态的耦合、角速度到欧拉角变化率的三角函数、以及空气阻力这类非线性项，形式很复杂。Koopman算子的核心想法是换个视角：不去线性化F本身，而是找一个观测函数g，让g(x_{k+1}) = K g(x_k)成立。这里的K是这个无穷维线性空间上的算子，它不近似F，而是精确地描述所有观测函数的演化。

这个「换空间」的代价是把有限维状态提升到无穷维函数空间。实际使用中我们用一组有限的字典函数{ψ_1, ψ_2, ..., ψ_N}张成一个子空间，在这个子空间里找一个K的有限维矩阵近似。只要字典选得好，非线性动力学在这个提升空间里可以被一个常数矩阵准确描述，且适用范围比在原始状态上的局部线性化宽得多。这正是Koopman算子跟传统泰勒展开线性化的本质区别：它不做局部近似，而是做全局嵌入。

实现上，我们把原始状态x提升为z = [x; sin(x_i); cos(x_j); x_i*x_j; ...]，然后假设z_{k+1} = A z_k + B u_k。这里的A矩阵就是Koopman矩阵的有限维实现，B矩阵处理控制输入。Koopman算子在这里起的作用不是某个具体数值，而是一整套提法：先定观测空间，再在观测空间上求线性演化。

2.2 EDMD算法：从快照数据回归Koopman矩阵

EDMD（Extended Dynamic Mode Decomposition，扩展动态模态分解）是求这个Koopman矩阵最常用的数据驱动方法。它的输入是一组快照对：状态序列x_1, x_2, ..., x_{m+1}，以及对应的下一时刻状态y_k = x_{k+1}。把每个快照都通过字典函数提升成向量ψ(x)，然后构造两块数据矩阵：

Ψ_X = [ψ(x_1), ψ(x_2), ..., ψ(x_m)]
Ψ_Y = [ψ(y_1), ψ(y_2), ..., ψ(y_m)]

如果字典函数张成的空间足够好，那么存在矩阵K使得Ψ_Y ≈ K Ψ_X。用最小二乘求解，K = Ψ_Y · pinv(Ψ_X)。等价写法是K = (G + λI)^{-1} A，其中G = Ψ_X Ψ_X^T是协方差矩阵，A = Ψ_X Ψ_Y^T是互协方差矩阵，λ是正则化系数。EDMD算法落地步骤非常直白：

1. 从仿真或实飞记录里取等间隔的状态序列与控制序列；
2. 定义字典函数集，写出提升函数lift(x)；
3. 对每个时刻k构造快照对(x_k, x_{k+1})，生成矩阵Ψ_X和Ψ_Y；
4. 用伪逆或带正则化的最小二乘求K；
5. 验证K的特征值是否位于单位圆内，以及提升模型的预测误差。

这个流程用Matlab实现，核心代码往往不超过二十行。但EDMD的结果质量高度依赖第2步的字典设计，字典太小无法表达非线性结构，字典太大会让回归矩阵病态。悬停点附近用几次多项式就够了，要覆盖大角度机动就必须引入三角函数与耦合项。

2.3 为什么偏偏是四旋翼：耦合、欠驱动与弱非线性假设

四旋翼是典型的欠驱动系统：四个电机输入控制六个自由度。它的动力学非线性主要体现在三个方面：推力矢量随姿态变化、欧拉角运动学方程中的三角函数耦合、以及高速飞行时的气动阻力。很多控制方案在小角度假设下线性化，把耦合项当作扰动交给PID去扛。这个方法在悬停和小机动场景够用，但做大速度前飞、大角度拉起时，线性模型里的耦合项误差会被MPC这类预测控制器放大。

数据驱动控制解决的是「模型精度与建模成本」的矛盾。Koopman算子对四旋翼特别合适的原因在于：四旋翼动力学不是强非线性或混沌系统，它的非线性主要来自姿态三角函数的耦合与二次型阻力，这些恰好都能被多项式字典和三角函数字典有效表达。也就是说，用EDMD把12维状态提升到30到80维的字典空间，线性模型的预测精度就足够支撑MPC设计。

另外一个工程层面的理由：传统方法需要精确的转动惯量、电机时间常数、升力系数等参数，这些参数在不同载荷条件下都会漂移。EDMD不需要这些物理参数的准确值，它直接从飞行数据里学出输入到状态演化的映射关系，换一架负载不同的机架时，重新录一段数据重跑一遍辨识就能更新模型。这个特性让它在快速原型验证里特别省事。

3. Matlab跑通EDMD辨识：从仿真数据到Koopman矩阵

3.1 在Matlab里搭一个带激励的四旋翼数据源

在Matlab里定义微分方程有几种常见姿势，可以用ode45配函数句柄，也可以用离散递推。做EDMD我一般直接用离散递推，因为后面辨识需要等间隔采样，递推模型天然就是一步预测的形式。下面是我在脚本里用的简化四旋翼刚体模型：速度环、位置环不展开，状态取为[位置, 速度, 欧拉角, 角速度]共12维。

% 四旋翼刚体动力学离散递推，dt为采样周期 function x_next = quad_dynamics(x, u, dt) % x: [px; py; pz; vx; vy; vz; phi; theta; psi; p; q; r] % u: [thrust; tau_phi; tau_theta; tau_psi] g = 9.81; m = 1.2; Ixx = 0.008; Iyy = 0.008; Izz = 0.014; phi = x(7); theta = x(8); psi = x(9); p = x(10); q = x(11); r = x(12); % 位置导数：机体速度旋转到世界系 R = euler_to_R(phi, theta, psi); vel_world = R * x(4:6); % 力方程：thrust沿机体z轴，重力沿世界z轴 accel_body = [0; 0; u(1)/m] - R' * [0; 0; g]; accel_world = R * accel_body; % 姿态运动学：欧拉角变化率与体角速度关系 phi_dot = p + q*sin(phi)*tan(theta) + r*cos(phi)*tan(theta); theta_dot = q*cos(phi) - r*sin(phi); psi_dot = (q*sin(phi) + r*cos(phi)) / cos(theta); % 角速度动力学：欧拉方程(忽略陀螺力矩) p_dot = (u(2) + (Iyy-Izz)*q*r) / Ixx; q_dot = (u(3) + (Izz-Ixx)*p*r) / Iyy; r_dot = (u(4) + (Ixx-Iyy)*p*q) / Izz; x_next = x + dt * [vel_world; accel_world; phi_dot; theta_dot; psi_dot; p_dot; q_dot; r_dot]; end

这段代码的R = euler_to_R(phi, theta, psi)负责把机体速度转到世界系，实际使用时可以用angle2dcm函数替代。关键参数是采样周期dt，EDMD要求数据是等间隔的，我习惯设dt = 0.01秒，也就是100Hz控制频率。状态序列x的维度是12，控制输入u的维度是4。升力系数、转动惯量在仿真里是「真实」参数，但EDMD辨识时完全用不到这些数值，这正是数据驱动的意义。

采集数据时不能只录悬停，激励信号很关键。我一般用悬停油门叠加扫频信号，频率从0.1Hz扫到10Hz，覆盖姿态和速度的主要动态范围。Matlab里生成这种激励用chirp函数即可，幅值取悬停油门值的5%到15%，太大容易让欧拉角超出小角度假设，太小又激发不出非线性特征。

3.2 字典函数怎么选：从单项式到三角函数再到耦合项

字典设计是EDMD里最影响结果的操作，也是我觉得最像「调参玄学」的地方。我的默认策略是：原始状态必须包含在字典里，否则最后预测出的状态没法直接从提升空间里读出来。除此之外，加三类扩展项——状态的平方项、欧拉角的三角函数项、以及速度与姿态角的耦合项。下面这个字典在12维状态上扩展到了54维，基本覆盖四旋翼在中等机动范围内的非线性特征。

% 构造提升向量z = lift(x)，字典由多项式、三角函数与耦合项组成 function z = lift(x) % 基础状态：位置、速度、欧拉角、角速度 z = x; % 速度平方项：表达阻力与动压效应 v = x(4:6); sq = v.^2; z = [z; sq]; % 姿态角三角函数：表达重力分量与旋转矩阵元素 sphi = sin(x(7)); cphi = cos(x(7)); stheta = sin(x(8)); ctheta = cos(x(8)); spsi = sin(x(9)); cpsi = cos(x(9)); z = [z; sphi; cphi; stheta; ctheta; spsi; cpsi]; % 速度与姿态角的耦合项：表达推力矢量在速度方向的投影 z = [z; v(1)*stheta; v(2)*sphi; v(3)*stheta; v(1)*cphi]; % 角速度平方项：表达陀螺进动力矩 omega = x(10:12); z = [z; omega(1)*omega(2); omega(2)*omega(3); omega(3)*omega(1)]; end

字典把12维状态提升到了54维。第一项到第十二项就是原始状态本身，这样预测结束后可以直接取z的前12个分量恢复出x的预测值。速度平方项用来拟合空气阻力；三角函数项覆盖欧拉角到旋转矩阵的映射；耦合项解决水平速度与姿态角之间的联动。这四个耦合项的选择依据是四旋翼的平移动力学：前飞时推力倾斜产生的水平加速度与theta成正比，侧向力与phi相关。

做EDMD最忌讳字典设计不看物理背景。比如把状态的所有五次方项都加进去，维数直接飙升到上百，条件数恶化不说，辨识出来的模型抗噪声能力也差。可以先用这个54维字典跑一遍，看预测误差主要出现在哪些状态上，再有针对性地补字典，而不是一次性堆满。

3.3 用伪逆求解Koopman矩阵：一条命令背后的数值风险

数据生成和字典定义好之后，EDMD求解本身只有几行代码。将状态序列提升成矩阵Ψ_X和Ψ_Y，然后用右除算子求解K。这个右除本质上是最小二乘，等价于pinv伪逆，但对矩阵结构更友好，速度也更快。

% 从仿真数据辨识Koopman矩阵 % u_all: 4xN 控制输入，x_all: 12xN 状态序列 N = size(x_all, 2); % 提升所有快照，生成字典矩阵 Psi_X = zeros(dim_z, N-1); Psi_Y = zeros(dim_z, N-1); for k = 1:N-1 Psi_X(:, k) = lift(x_all(:, k)); Psi_Y(:, k) = lift(x_all(:, k+1)); end % 最小二乘求Koopman矩阵，右除等价于 Psi_Y * pinv(Psi_X) K = Psi_Y / Psi_X; % 检查K的特征值是否都在单位圆内 ev = eig(K); disp(['max abs eigenvalue = ', num2str(max(abs(ev)))]);

这段代码里Psi_X的每一列是一个快照的提升向量，Psi_Y的对应列是下一时刻的提升向量。K的维度是dim_z乘dim_z，这里就是54乘54。如果max abs eigenvalue大于1.01，说明模型有不稳定模态，通常是字典病态或数据激励不足导致的，要先排查而不是急着接控制器。

单靠特征值还不够，我会再做一次开环预测验证：从某个初始状态出发，用辨识出的K独立推50步，把预测轨迹与真实仿真轨迹对比，计算归一化均方根误差。

% 开环预测验证：对比EDMD模型与真实模型的轨迹 x_pred = zeros(12, H); z = lift(x_all(:, 1)); % 取第一帧作为初始提升状态 for k = 1:H z = K * z; x_pred(:, k) = z(1:12); % 前12维还原原始状态 end % 归一化误差 err = x_pred - x_all(:, 1:H); nmse = sum(err.^2, 'all') / sum(x_all(:, 1:H).^2, 'all'); fprintf('NMSE = %.4f\n', nmse);

预测误差如果超过5%，我一般先怀疑字典不充分，比如缺少某个耦合项；如果误差集中在姿态角，优先在字典里补三角函数相关项；如果误差在所有维度均匀偏大，则考虑数据激励幅值不够。开环验证通过后，再进入控制器设计阶段。

4. 在Koopman线性模型上设计MPC：预测控制与Matlab求解

4.1 从提升状态到线性预测模型：A、B矩阵的组合回归

上一章的K只表达了无控制输入的自由演化。四旋翼必须有控制输入，所以需要把EDMD扩展成带输入的回归形式：z_{k+1} = A z_k + B u_k。严格来说，四旋翼是控制仿射系统，提升之后理论上会变成双线性结构（输入与状态相乘的项），但工程实践里我们在工作点附近把B当作常数矩阵处理，悬停和中低速机动范围内精度足够。

求A和B的做法是把输入也放进回归变量里：一边是[Ψ_X; u]，另一边是Ψ_Y，这样解出来的矩阵前半部分对应A，后半部分对应B。

% 带控制输入的EDMD：z_{k+1} = A * z_k + B * u_k % 构造组合回归矩阵 X_reg = [Psi_X; u_all(:, 1:N-1)]; % (dim_z + 4) x (N-1) Y_reg = Psi_Y; % dim_z x (N-1) % 最小二乘组合求解 M = Y_reg / X_reg; A = M(:, 1:dim_z); B = M(:, dim_z+1:end); % 稳定性检查 if max(abs(eig(A))) > 1.01 warning('A矩阵存在不稳定模态，请检查字典与激励信号'); end

这里B矩阵的物理含义是：每个控制通道（推力、滚转力矩、俯仰力矩、偏航力矩）对全部54个字典状态的单步影响。因为B是从数据里回归出来的，它天然包含了输入到加速度、角加速度的增益关系，不需要单独辨识转动惯量。要注意的是B是常数矩阵，一旦机动幅度变大，B的精度会下降，这是这种方法的边界。

4.2 用quadprog写线性MPC：约束、参考值与代价函数

有了线性模型z_{k+1} = A z_k + B u_k，MPC设计就成了标准的线性二次规划问题。目标函数惩罚三类量：未来预测状态与参考值的偏差、控制量大小、以及控制量变化率。约束条件包括电机推力上下限和姿态角速率限制。在Matlab里我用quadprog求解这个QP问题。

function [u_opt, z_pred] = koopman_mpc(z0, u_prev, x_ref, A, B, C, Np, Nc, Q, R, umax) % z0: 当前提升状态，C: 从提升状态提取原始状态的矩阵 % Np: 预测时域步数，Nc: 控制时域步数 dim_z = length(z0); dim_u = size(B, 2); % 构造预测矩阵：Z_pred = Phi_z * z0 + Phi_u * U % 这里用循环构造，方便理解；大规模场景可以换为稀疏矩阵 Phi_z = zeros(Np*dim_z, dim_z); Phi_u = zeros(Np*dim_z, Nc*dim_u); % 状态传播：z_{k+1} = A z_k + B u_k A_pow = eye(dim_z); for i = 1:Np A_pow = A * A_pow; Phi_z((i-1)*dim_z+1 : i*dim_z, :) = A_pow; for j = 1:min(i, Nc) A_pow_ij = A^(i-j); Phi_u((i-1)*dim_z+1 : i*dim_z, (j-1)*dim_u+1 : j*dim_u) = A_pow_ij * B; end end % 目标函数矩阵：只惩罚原始状态分量(C*z)，不需要给所有字典状态设参考 Q_bar = kron(eye(Np), C' * Q * C); R_bar = kron(eye(Nc), R); H = Phi_u' * Q_bar * Phi_u + R_bar; f = Phi_u' * Q_bar * (Phi_z * z0 - repmat(x_ref, Np, 1)); % 输入幅值约束 lb = -umax * ones(Nc*dim_u, 1); ub = umax * ones(Nc*dim_u, 1); % 求解QP options = optimoptions('quadprog', 'Display', 'off'); U = quadprog(2*H, 2*f, [], [], [], [], lb, ub, [], options); u_opt = U(1:dim_u); z_pred = Phi_z * z0 + Phi_u * U; end

这个MPC函数里有几个容易绕进去的点。第一个是参考值只能针对原始状态定义，比如位置、速度、姿态角，不能针对全部54维字典状态，因为RBF或平方项根本没有物理含义。所以我用C矩阵把预测状态映射回12维原始状态，再用Q矩阵加权。第二个是构造Phi_u矩阵时j循环只跑到min(i, Nc)，控制时域之后的输入保持为0，这是MPC约定俗成的处理。第三个是quadprog被写成了2H和2f，这是因为quadprog默认最小化0.5*x'Hx + f'x，而标准MPC代价是x'Hx + f'x，要换算系数。

4.3 闭环仿真与控制器参数调法

闭环仿真就是把MPC编进主循环里，每步采样后用当前状态求解控制量，再作用到四旋翼模型上。这是在Matlab里验证整套方案的最后一步，也是暴露模型失配问题的关键时刻。

% 闭环仿真主循环 x = x_init; u = u_hover; x_log = zeros(12, T_total); for k = 1:T_total z_k = lift(x); [u, ~] = koopman_mpc(z_k, u, x_ref, A, B, C, 30, 5, ... diag([10 10 10 1 1 1 5 5 5 0.1 0.1 0.1]), ... 0.1*eye(4), 2.0); x = quad_dynamics(x, u, dt); x_log(:, k) = x; end

预测时域Np我习惯取30步，对应0.3秒的预测长度，对四旋翼姿态控制足够，又不至于让QP矩阵过大。控制时域Nc取5到8步，太大会让优化变量过多、求解变慢，太小会丧失MPC的预见性。Q权重里位置和速度的权重是关键：想让无人机快速跟踪轨迹就加大位置权重，想让姿态过渡平顺就提高姿态角权重并压低角速度权重。R取0.1 * eye(4)表示对四个控制通道施加相同程度的惩罚，实际调试如果发现油门抖动，先把R对应的推力通道权重加大到0.5。

调参顺序上我的经验是：先固定R，把Q从大到小试一遍，找到临界稳定的点；再固定Q，微调R消除高频抖动；最后再微调Np和Nc。如果闭环一开始就发散，先别急着调权重，回头检查A矩阵特征值和开环预测误差，大概率是模型本身没辨识好。

5. EDMD落地最容易翻车的5个坑：从字典病态到闭环发散

5.1 数值病态与字典设计问题：条件数过大导致K矩阵不可用

现象：Psi_X矩阵的条件数高达1e14，求解出的K不满足稳定性检查，特征值散落在单位圆外，开环预测几步就爆掉。

原因：字典里混入了量级差异极大的特征。比如位置状态可能是10米量级，角速度可能是2弧度每秒量级，而速度平方项直接到了100。这些不同量纲的观测函数放在同一个矩阵里，最小二乘会被大数值项主导，小数值项对应的动力学被数值噪声淹没。另一个常见诱因是两个字典函数高度线性相关，比如同时加入了sin(phi)和phi，或者同时加入了cos(theta)和1 - theta^2/2。

解决：先对每个字典分量做归一化。统计训练数据里每个字典分量的均值和标准差，把Psi_X和Psi_Y都减去均值再除以标准差，求解完成后再把A矩阵还原到原始量纲。用Matlab的zscore函数就能做。同时检查字典函数是否线性独立，用rank(Psi_X)对比矩阵行数，缺秩就删掉冗余项。我常用的做法是先用54维基础字典，出现问题再逐步加项，而不是一次性堆到上百维。

5.2 欧拉角在正负180度跳变导致预测发散

现象：开环验证里俯仰角或偏航角预测误差随时间线性增大，一旦无人机做了超过90度的翻转，预测立刻跳到完全相反的方向。

原因：欧拉角本身不是连续量，psi从+179度变到-179度时数值跳了358度，但实际飞机姿态只变了2度。EDMD回归的是数值映射关系，它无法理解这种跳变背后的拓扑含义，于是把一次正常的跨越边界当作一次剧烈的动力学事件，预测自然爆掉。

解决：不要直接用欧拉角作为字典分量，改用旋转矩阵的元素或四元数。最省事的方法是把姿态表示换成欧拉角的三角函数组合，也就是字典里的sphi、ctheta这些项，原始状态里的欧拉角仍然保留给参考值和显示用，但预测循环内部操作的是提升状态。如果一定要直接预测欧拉角，就在采集数据之前把psi序列做unwrap处理，让角度序列连续变化，但这个方法只对小角度范围有效。

5.3 数据激励不足：悬停数据辨识出的模型在机动时失效

现象：用悬停加小幅扰动数据辨识出的模型，开环验证NMSE只有1%，但接上MPC做大机动跟踪时，预测轨迹和实际轨迹迅速分离，控制器在机动指令发出后出现明显的滞后和振荡。

原因：EDMD本质上是在拟合字典空间内的输入输出映射。训练数据没有覆盖的动力学区域，模型完全靠外推。悬停数据里速度始终接近零，速度平方项的系数拟合出来接近于随机值；一旦做大速度前飞，平方项被激活，错误的系数就直接体现在预测误差里。这跟神经网络过拟合是一个道理，只是EDMD的字典是显式的，更容易看清症结。

解决：采集数据时把激励幅值放宽到目标工作范围。做轨迹跟踪就往对应的速度区间扫频，做大姿态机动就加入方波形式的姿态指令激励。更稳妥的做法是设计一个递进式的数据采集流程：先用小幅激励辨识一个粗模型，在这个模型上跑一个保守的MPC，再用闭环数据二次辨识，逐步扩大工作范围。用闭环数据做EDMD时要注意把输入也记录下来，闭环EDMD的关键是回归新数据时把历史模型作为先验正则项，防止新数据覆盖掉悬停附近的动态。

5.4 闭环里预测误差累积：开环验证通过但闭环发散

现象：开环预测50步的NMSE只有3%，看起来模型质量很好。一旦接入MPC闭环，第一个步长还正常，二十步之后控制量开始高频振荡，最终发散。

原因：开环预测误差是单步误差不断累加的结果，而闭环控制会实时纠正真实状态与预测状态的偏差。问题往往不在模型精度，而在于MPC对模型失配的反馈增益太高。Koopman模型在字典子空间内的预测是线性的，但真实系统的非线性残留会让MPC的优化解在某些状态上产生过大的控制增量，这些增量反过来又让模型在下一步更不准确，形成恶性循环。

解决：给MPC增加鲁棒性手段。第一是在代价函数里加入控制增量惩罚项，限制控制量每一步的变化幅度，避免模型失配带来的控制量突变。第二是把持续扰动当作附加状态：在预测模型里加一项L(z - z_true)作为误差反馈修正项，其中L通过对比历史预测误差统计得出。第三是保守地缩短预测时域Np，从30步收到15步，虽然预见性变弱，但模型失配的影响在短时域内被摊薄。最后一种做法是牺牲一部分性能把Q权重里速度项降低，让控制器不那么激进地响应预测偏差。

5.5 求解方式用错：右除、伪逆与截断SVD的选择

现象：同样的快照数据，用Psi_Y / Psi_X快速求解得到的结果正常，换成pinv(Psi_X) * Psi_Y之后，K的特征值出现变化，MPC效果也差了一截。

原因：右除和pinv在线性代数里数学上等价，但数值行为不同。右除算子对欠定和接近病态的矩阵会走最小二乘路径，遇到奇异值极小的方向会自动丢弃部分信息。pinv默认保留所有非零奇异值，当Psi_X存在接近零的奇异值时，逆运算会把噪声放大到离谱的程度。两种求解方式在不病态的矩阵上结果一致，但四旋翼字典矩阵经常不是良态的。

解决：统一使用带截断的求解方式。先对Psi_X做奇异值分解，把小于最大奇异值千分之一的奇异值直接归零，再求伪逆。Matlab里可以用[U,S,V] = svd(Psi_X)，把S的对角元小于eps * max_sv的部分置零，然后重建伪逆。或者更省事的方式是给正规方程加一个小正则项。这一步的差别在开环验证里不容易看出来，但往往决定了闭环能不能稳定。

6. 从离线辨识到在线修正：收敛性检查与字典再设计

Koopman模型不是辨识完就能一直用的。每次换载荷、改机架、换旋翼，动力学都会变，数据驱动模型最实际的维护方式是新数据进来之后做增量更新。EDMD的增量版实现起来很简单：EDMD求解的中间量G和A可以累加，新数据到来时把新的快照对加入G和A，不需要重新遍历历史数据。Matlab里维护这两个累加矩阵，用G_new = G_old + Psi_X_new * Psi_X_new'，A_new = A_old + Psi_X_new * Psi_Y_new'，再对G_new加一个极小的单位阵正则项，直接解K = G_new \ A_new。

字典再设计要基于残差分析而不是拍脑袋。完成一轮辨识后，把预测残差按12个原始状态分量分别统计，哪个方向误差大就在哪个方向补字典。比如位置预测误差偏大，通常是速度与姿态的耦合项不足；偏航角误差偏大，通常是psi的三角函数项缺失或者欧拉角跨越边界。每次补字典之后重新跑开环预测，观察NMSE是否真的下降，如果加了新项但误差没变，说明该项与现有字典线性相关，直接去掉，不要舍不得。

我这几年的习惯是：每个项目都固定保留两份脚本，一份是EDMD辨识与验证，一份是MPC闭环仿真。辨识脚本里把字典结构、激励信号、求解方式全部参数化，换新机型时只改参数不重写代码。闭环仿真里固定记录特征值、开环NMSE和闭环跟踪误差三个指标，任何修改都以这三个指标不变差为准。先开环后闭环，先调字典后调控制器，这个顺序帮我避开了大半的调试翻车现场。这套方案从仿真到半实物仿真的迁移路线也清晰：把辨识脚本里的仿真数据源替换成飞行日志，MPC的输入输出接口不变，就可以直接验证真实飞行数据下的控制效果。希望帮到你。

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