做控制方向的朋友,大概率都经历过这种尴尬:你搭好一个非线性模型预测控制器,离线仿真效果漂亮,一上实时环境,每一步要解的非线性非凸优化动辄几十毫秒,采样周期根本撑不住。我之前在一个倒立摆Demo上就卡在这儿——系统只是简单的两状态模型,但非线性项加上约束之后,在线求解时间始终压不下来,控制器参数改来改去也只能勉强跑在“能跑但不稳”的边缘。
后来我把思路换成了Koopman模型预测MPC,用数据驱动的方式先把这个非线性受控系统“提升”到一个高维线性空间,得到一个线性预测器,再把MPC问题退化成标准的二次规划。同样的问题,quadprog一跑就是几毫秒,精度反而比原来的非凸求解更高。这篇文章就围绕这个方案展开:讲清楚Koopman算子为什么能把非线性系统变成线性预测器、如何用EDMD离线辨识模型,再给出完整的Matlab代码实现,最后聊聊我实测中踩过的坑。
这套内容适合正在做机器人控制、过程控制,或者刚接触Koopman理论、想在Matlab里快速验证MPC效果的人。不需要太高深的数学基础,但至少要知道什么是状态方程、什么是二次规划,并且愿意亲手跑代码。我尽量把“为什么这么做”和“到底怎么做”两条线都讲透。
1. 非线性MPC的困局:为什么我想把系统“变线性”
1.1 MPC的本质与非线性痛点
模型预测控制的核心并不复杂:在当前时刻拿到状态x(k),用系统模型向前预测一段时域内的状态轨迹,然后在线求解一个带约束的最优控制问题,得到未来几步的控制序列,只把第一步作用到系统上,下一时刻再滚动重复。它的吸引力在于显式处理约束,非常适合工业过程、机器人运动控制这类对安全性要求高的场景。
问题出在预测模型上。如果系统是线性的,比如x_{k+1}=Ax_k+Bu_k,那么整个预测是线性的,代价函数是二次的,约束是线性的,优化问题就是一个标准二次规划,Matlab里quadprog、OSQP这类求解器可以非常快地解完。但如果系统本身是非线性的,比如倒立摆里的sin项、机械臂里的科式力和离心力项,为了得到精确的预测轨迹,每一步都要在线求解一个非线性非凸优化。
非线性非凸优化有多贵?贵在它需要反复迭代求解KKT条件,用内点法或序列二次规划(SQP)都需要多次梯度计算和黑塞矩阵更新。对于一个小规模倒立摆,NLP求解普遍需要几十毫秒甚至更久,而很多控制系统的采样周期只有5到20毫秒。这意味着单纯堆算力很难解决实时性问题,必须从模型结构上做文章。
1.2 传统线性化路线的三条“折中方案”及其局限
既然非线性MPC在线求解太贵,工程上最常见的思路就是“想办法把模型变线性”。我实测过三类折中方案,各有各的局限。
第一类是在工作点附近做雅可比线性化。系统在工作点x0附近被写成小偏差线性模型,然后对这个线性模型做MPC。这种方案在小偏差工况下效果很好,但问题很拧巴:每换一个工作点就要重新计算一次雅可比矩阵、重新预测矩阵,而且一旦状态偏离工作点较远,模型误差急剧增大。控制器会“近视”,看不了太远的未来,预测时域稍微拉长一点就会出现预测轨迹漂移。
第二类是反馈线性化,通过状态变换和一个精心设计的控制律,把非线性系统中的非线性项直接抵消掉,得到一个等效的线性系统。理论上非常漂亮,但需要系统满足相对阶等严格条件,必须精确知道系统模型的全部非线性结构,稍微有一点模型不确定性,闭环性能就大打折扣。
第三类是显式MPC,把MPC的控制律预先离线分区域求出,在线退化成查表。这在小规模问题里很实用,但状态空间会被划分成大量多面体区域,状态维度稍微一高、约束稍微一多,离线运算量就爆炸式增长。
这三条路线本质上都在做同一件事:牺牲非线性特性,换取在线计算效率。差别只是用哪种方式“近似”。我的体会是,它们各自的牺牲代价太明显——要么只在小范围内有效,要么对模型精度要求苛刻,要么离线成本高。我一直希望找到一个“更体面”的线性化方式。
1.3 Koopman MPC的破局思路:换一个空间看系统
Koopman MPC让人眼前一亮的地方在于,它不直接对非线性状态方程做线性化,而是切换到一个更高维的“观测函数空间”去描述系统演化。
通俗地说:原来你在原始状态空间里用x_{k+1}=f(x_k,u_k)描述系统,这个f可能是弯弯绕绕的非线性映射。Koopman的思路是,如果你构造一组新坐标(观测函数),比如ψ(x)=[x; sin(x); cos(x); x²],让系统在这组新坐标下的演化变成线性的,即ψ(x_{k+1})≈Aψ(x_k)+Bu_k,那么在这个提升后的高维空间里,非线性预测问题就重新变成了线性预测问题,MPC也重新回到了线性二次规划的舒适区。
这个思路听起来有点像“换坐标系看问题”。就像同一个物理过程,在直角坐标系里是一条曲线,换到合适的参数坐标系里可能就变成了一条直线。Koopman理论给这件事提供了数学基础:只要能找到足够丰富的观测函数集合,非线性系统的演化在观测函数空间里可以用线性算子描述,代价只是算子维度可能无限大。实际工程中,我们截取一个有限维逼近,这就是整个方法的核心。
2. Koopman算子到底在干什么:从观测函数到高维线性空间
2.1 一个公式看懂Koopman算子
先看一个不带动力的离散系统x_{k+1}=f(x_k)。传统方式关心状态本身如何演化。Koopman算子的关注点变了,它关心的是“任意的观测函数g(x)”如何演化:
(Kg)(x) = g(f(x))
这里的g可以是一个标量函数,比如“摆角的正弦值”g(x)=sin(x₁),也可以是向量值函数。Koopman算子K的作用,就是把“当前状态的观测值g(x)”映射到“下一个状态的观测值g(f(x))”,也就是一步时间演化在观测函数空间里的体现。
关键在于,K本身是一个线性算子,哪怕底层的f非常非线性。因为算子的线性体现在对函数进行线性组合时满足L(g₁+g₂)=Lg₁+Lg₂,这一步直接可以验证。这给了我们一个非常重要的结论:非线性系统的状态演化和观测函数演化之间,存在一层“线性结构”,只是它生活在无穷维的函数空间里,而不是原始的有限维状态空间。
2.2 受控系统与有限维近似形式
对受控系统x_{k+1}=f(x_k,u_k),Koopman理论也可以推广,但工程实现上更常见的是直接假设存在一组字典函数ψ(x),使得:
ψ(x_{k+1}) ≈ Aψ(x_k) + Bu_k
也就是把控制输入作为外生线性项加进来,把状态提升空间里的递推关系近似成线性状态空间模型。这里的A是一个n_lift×n_lift的矩阵,B是n_lift×n_u的矩阵,n_lift是字典函数个数,也就是提升空间的维度。
严格来说,这只是对受控Koopman的一种简化近似,并不能覆盖所有非线性受控系统。但在工程实践中,如果字典函数选得足够好、训练数据覆盖得当,这种形式已经足够支撑MPC的预测需求。在做EDMD辨识时,我们实际上就是在用数据拟合A和B,而不是从系统的显式方程去推导Koopman算子。
2.3 为什么“无限维线性”能截断成“有限维近似线性”
真正让我第一次读懂Koopman理论的,是“坐标基底”这个类比。原始状态就好比你在看一个球体的三维投影,你看到的是一个圆在动,轨迹是弯曲的。Koopman算子把问题放到无穷维的函数空间里,相当于加了很多很多坐标轴,每个坐标轴都从不同角度观测系统状态。在这些坐标轴张成的空间里,系统的运动方程反而变得简单。
但无穷维没法直接算,所以工程上必须做截断:只选取n_lift个字典函数,把观测函数空间限制在由这些基函数张成的有限维子空间上。Koopman算子对空间的演化,在这个子空间上的投影就是矩阵A。截断必然带来近似误差,但有一个好消息:如果字典函数取得足够“贴近”系统的动力学特征,误差可以被压到很小。经典的Koopman模态理论告诉我们,系统响应中最重要的几个模态,通常可以对应低频的主要特征,用几十个字典函数往往就能描述得相当好。
我自己的体会是,有限维截断这件事才是整个方法的核心艺术。因为数学上“存在某个无穷维线性算子”是优雅的理论,但在Matlab里跑通MPC,你需要的是一个能用的矩阵A。字典越多,逼近越好,但维度太高又会让EDMD出现过拟合和数值不稳定。这个尺度怎么拿捏,我后面专门用一节来讲。
3. 把理论变成模型:EDMD数据驱动的辨识流程
3.1 EDMD到底在做什么
EDMD(扩展动态模态分解)本质上是数据驱动地找一个矩阵A和B,使得提升空间里的线性递推关系能够最好地拟合从数据中观测到的非线性系统演化。给定一批数据对(x_k, x_{k+1}, u_k),我们希望最小化:
min ||ψ(x_{k+1}) - (Aψ(x_k) + Bu_k)||²
把所有样本堆起来,就是一个线性最小二乘问题,可以一次性求出A和B。这就是EDMD的全部核心,没有更多花哨的东西。让我直白地说:它就是把非线性系统在提升空间里“线性回归”了一下。
不过这个“线性回归”有一个前提,就是要设计好特征(即字典函数ψ),还要保证数据覆盖范围足够宽。特征设计得好不好,直接决定回归出来的A矩阵能不能真正捕捉到非线性特性;数据覆盖够不够,直接决定回归结果能不能泛化到目标工作域。
3.2 字典函数怎么选:倒立摆案例分析
字典函数的选择没有放之四海而皆准的公式,但有几条可操作的“套路”。我在倒立摆案例里采用的字典,基本上是“原始状态 + 非线性项 + 交叉项”的组合:
- 原始状态x₁、x₂,保证提升空间中前两个维度能直接映射回原始状态;
- 三角函数项sin(x₁)、cos(x₁)-1,用来捕捉摆角的大范围非线性特性;
- 更高阶项sin(2x₁)、cos(2x₁)-1,补充更高频率的模态信息;
- 二次多项式项x₁²、x₂²、x₁x₂,用来表达状态之间的耦合和局部非线性。
有一点值得强调:字典函数中务必保留原始状态本身作为前几个维度,这样后续MPC做约束时,可以把约束直接施加在提升空间的前几个分量上,轻松映射回现实世界。如果字典完全由抽象特征构成,约束处理会非常别扭。
对于其他系统,RBF核函数和随机傅里叶特征也是常见选择。RBF的好处是通用,不需要太多先验知识,但维度稍微高一点就容易病态;傅里叶特征适合处理强振荡行为。我个人的偏好是先从物理直觉出发,把明显非线性模态加进字典,如果效果不好,再用数据驱动方式(比如Koopman特征函数逼近)去补充。
3.3 岭回归求解与验证集检查
矩阵拟合时,直接最小二乘有风险。因为提升空间维度一旦超过几十维,数据矩阵很容易出现共线性,导致伪逆结果数值不稳。我的做法是加一个小的正则项,用岭回归来解:
K = Ψ_y Θ^T (ΘΘ^T + λI)^{-1}
这里的Θ是提升状态和控制输入堆叠成的矩阵,Ψ_y是下一时刻的提升状态矩阵。λ通常在1e-6到1e-3之间,具体值需要通过验证集来确定,不宜拍脑袋随便设。
训练完之后,千万别看个训练误差就说模型好。我每次都会专门留出一段“验证数据”,单独算一步预测误差和滚动预测误差。如果训练数据上的拟合误差很小,但验证数据上误差很大,说明过拟合了,需要降低字典规模或者加大正则项。如果两者误差都大,多半是字典函数不够丰富,或者训练数据激励不充分。
4. Matlab代码实现:从倒立摆仿真到线性MPC闭环
这一节我直接给可运行的Matlab代码脉络,方便你在自己的机器上复现。我用一个受控倒立摆作为测试对象,状态是摆角x₁和摆角速度x₂,控制输入是电机力矩u。离散时间模型用欧拉近似,采样周期取dt=0.02秒。
%% 倒立摆系统定义(离散时间欧拉近似) g = 9.81; L = 0.5; m = 0.3; b = 0.05; dt = 0.02; nx = 2; nu = 1; f = @(x, u) [... x(1) + dt*x(2); ... x(2) + dt*(g/L*sin(x(1)) - b/(m*L^2)*x(2) + u/(m*L^2)) ... ];4.1 训练数据生成:激励信号要覆盖目标工作域
训练数据的质量直接决定了Koopman模型的好坏。我用带限幅的随机方波作为激励信号,方波在每个幅值上保持若干个采样周期,这样比纯高斯白噪声更容易激发系统的大范围动态特性,也更接近真实控制信号的形态。
%% 训练数据采集 rng(7); N = 40000; % 训练样本数 X = zeros(nx, N); % 状态历史 U = zeros(nu, N-1); % 控制输入历史 x = [0; 0]; u_amp = 1.5; % 输入幅值上限 hold_steps = 10; % 方波保持步数 for k = 1:N-1 if mod(k, hold_steps) == 1 u_k = u_amp * (2 * (rand > 0.5) - 1); end U(:, k) = u_k; X(:, k) = x; x = f(x, u_k); end需要注意,初始状态我取在了倒立摆的下垂平衡点附近(角度为0,角速度为0)。由于倒立摆模型的sin项在0附近仍是非线性的,训练数据覆盖±1.5力矩激励下的动态范围,已经足够让Koopman模型学会向上摆起和稳定的动力学。如果你想让MPC控制到竖直向上平衡点,最好在采集数据时把初始状态放在竖直向上附近,或者覆盖包括竖直向上邻域在内的整个工作范围。实践下来,训练数据范围覆盖目标工作点这一点再怎么强调都不过分。
4.2 构建字典函数并拟合Koopman矩阵
字典函数实现上有一个细节:为了保证批量计算时矩阵每一列都能通过向量操作完成,我直接用一个匿名函数来表示字典,并对输入做归一化和模式匹配。下面的lift_fun返回一个列向量,输入是2维状态。
%% 字典函数 lift_fun = @(x) [ 1; % 常数项 x(1); x(2); % 原始状态 sin(x(1)); cos(x(1))-1; sin(2*x(1)); cos(2*x(1))-1; x(1)^2; x(2)^2; x(1)*x(2) ]; n_lift = length(lift_fun([0;0]));这里用cos(x)-1主要是为了在x=0时让该字典项也为0,避免和常数项产生共线性。如果你直接放cos(x),也没错,但会让回归问题稍微复杂一点,因为cos(0)=1等于常数项取值。
拟合A和B的代码如下:
%% 构造提升矩阵 PSI_X = zeros(n_lift, N-1); PSI_Y = zeros(n_lift, N-1); for k = 1:N-1 PSI_X(:, k) = lift_fun(X(:, k)); PSI_Y(:, k) = lift_fun(X(:, k+1)); end Theta = [PSI_X; U]; lambda = 1e-5; Kmat = PSI_Y * Theta' / (Theta * Theta' + lambda * eye(n_lift + nu)); A_koop = Kmat(:, 1:n_lift); B_koop = Kmat(:, n_lift+1:end);这段代码的核心就一行:Kmat = PSI_Y * Theta' / (Theta * Theta' + λI)。它同时解决了A和B的拟合,因为Theta把提升状态和控制输入并排堆在一起,最小二乘解一次得到两个矩阵。注意Matlab里/是右除,等价于乘以伪逆。lambda加在分母矩阵的主对角线上,起到正则化作用,防止Theta的协方差矩阵不可逆。
4.3 构建MPC预测矩阵:把多步预测写成线性映射
MPC的核心是把“未来Np步的状态轨迹”显式表达成“当前提升状态z0”和“未来控制序列U”的线性函数:
Z = Φ·z0 + Θ·U
其中Z是Np*n_lift维的大向量,包含从第1步到第Np步的提升状态;Φ和Θ是预先根据A_koop、B_koop构建好的预测矩阵。这一步构建必须仔细,它是后续所有QP代码的基础。
%% 构建预测矩阵 Phi 和 Theta_mat function [Phi, Theta_mat] = build_prediction_matrices(A, B, Np) n_lift = size(A,1); n_u = size(B,2); Phi = zeros(Np*n_lift, n_lift); Theta_mat = zeros(Np*n_lift, Np*n_u); Atmp = eye(n_lift); for k = 1:Np Atmp = A * Atmp; Phi((k-1)*n_lift+1 : k*n_lift, :) = Atmp; for j = 0:k-1 A_pow = A^(k-1-j); Theta_mat((k-1)*n_lift+1 : k*n_lift, ... j*n_u+1 : (j+1)*n_u) = A_pow * B; end end end这个函数返回的Theta_mat就是预测矩阵,它的第k行块保存了第k步预测输入对最终预测状态的影响系数。因为A^(k-1-j)B本质上描述的是“在第j步施加的输入u(j)如何通过剩余的时间步影响第k步的状态”,这个因果结构用矩阵乘法一笔写就,非常清晰。
4.4 把MPC写成二次规划:quadprog就是求解器
有了预测矩阵之后,MPC变成了标准的QP问题。决策变量是未来Np步的控制序列U,目标函数是预测状态与参考状态的二次偏差加上控制能量惩罚:
min 0.5·U'·H·U + f'·U 约束:U_min ≤ U ≤ U_max,以及状态范围约束(通过提升空间的前nx个分量实现)
具体代码中,我先把权重矩阵扩展成预测时域上的分块对角阵,然后构建H、f、线性不等式约束Aineq/bineq。
%% MPC权重与约束设置 Np = 30; % 预测时域 Q = blkdiag(10, 1, 0.1*eye(n_lift-2)); % 只惩罚前两维(真实状态) R = 0.1 * eye(nu); Q_bar = kron(eye(Np), Q); R_bar = kron(eye(Np), R); [Phi, Theta_mat] = build_prediction_matrices(A_koop, B_koop, Np); H = Theta_mat' * Q_bar * Theta_mat + R_bar; H = (H + H') / 2; % 强制对称,避免数值误差破坏QP凸性 x0 = [0.6; 0.2]; % 初始摆角约35度 z0 = lift_fun(x0); z_ref = lift_fun([pi; 0]); % 目标:竖直向上平衡点(角度=pi) % 注意:实际中根据你定义的角度坐标系调整参考值,我这里的0是下垂平衡点 %% 计算线性项 f Z_ref_bar = kron(ones(Np,1), z_ref); f_vec = 2 * Theta_mat' * Q_bar * (Phi * z0 - Z_ref_bar); %% 状态约束:通过 C*z 将原状态拉回 C = [eye(nx), zeros(nx, n_lift-nx)]; C_bar = kron(eye(Np), C); x_min = -0.5*ones(nx,1); x_max = 0.5*ones(nx,1); % 此处是示意:倒立摆摆角约束和角速度约束可按工程需求设定 Aineq = [C_bar * Theta_mat; -C_bar * Theta_mat]; bineq = [kron(ones(Np,1), x_max) - C_bar*Phi*z0; -kron(ones(Np,1), x_min) + C_bar*Phi*z0]; %% 控制输入约束 umax = 1.5 * ones(nu,1); umin = -1.5 * ones(nu,1); lb = kron(ones(Np,1), umin); ub = kron(ones(Np,1), umax); %% 求解QP并取第一步控制 options = optimoptions('quadprog', 'Display', 'off', 'Algorithm', 'interior-point-convex'); [u_seq, ~, exitflag] = quadprog(H, f_vec, Aineq, bineq, [], [], lb, ub, [], options); if exitflag > 0 u_current = u_seq(1:nu); else warning('QP求解失败,使用上一步控制量'); u_current = 0; end这段代码最需要注意的有两点。第一,z_ref的设置必须和你定义的字典函数一致,如果字典里包含常数项1,那么z_ref里对应的常数项也必须是1,否则参考信号会引入一个永远无法消除的偏置。第二,状态约束只应该施加在原始状态对应的维度上,也就是通过C矩阵把提升状态的前两个分量映射回来,而不是直接约束全部n_lift个提升维度——那些抽象维度没有物理意义,约束它们只会让QP变得难以求解。
4.5 闭环主循环:滚动预测,每次只施加第一步
把上面的步骤包成一个函数,然后放进闭环循环里:
%% 闭环仿真 Tsim = 3.0; % 仿真时长 steps = round(Tsim / dt); x = [0.6; 0.2]; % 初始摆角约35度 u_log = zeros(nu, steps); x_log = zeros(nx, steps+1); x_log(:,1) = x; for k = 1:steps u_current = koopman_mpc_controller(x, z_ref, A_koop, B_koop, ... Q, R, Np, umax, x_min, x_max, lift_fun); u_log(:,k) = u_current; x = f(x, u_current); % 作用到真实非线性系统 x_log(:,k+1) = x; endkoopman_mpc_controller就是上一小节代码封装后的函数。实测下来,每个控制周期里quadprog求解QP的耗时在2到5毫秒左右,远低于采样周期20毫秒,实时性完全够用。
5. 实测效果与参数取舍:Koopman能走多远
5.1 与工作点线性化MPC的对比
我同时实现了两个版本:一个是在竖直向上平衡点附近做雅可比线性化的线性MPC,另一个是Koopman MPC。控制对象完全一样,都是同一个倒立摆仿真模型。对比情况如下表:
| 初始摆角偏差 | 工作点线性化MPC | Koopman MPC |
|---|---|---|
| 5度 | 稳定,调节时间约0.6秒 | 稳定,调节时间约0.5秒 |
| 15度 | 勉强稳定,有明显超调 | 稳定,调节时间约0.8秒 |
| 35度 | 发散,无法镇定 | 稳定,调节时间约1.5秒 |
| 60度 | 发散,预测完全失真 | 基本稳定,但靠近训练数据边界,有轻微振荡 |
这个结果直观地说明,工作点线性化模型的有效控制范围大约只在平衡点附近±10度内;超出这个范围后,模型预测失真导致控制器无法有效镇定。Koopman MPC因为是在提升空间里逼近了全局非线性动态,所以在大偏差工况下依然能保持预测精度,控制范围宽得多。
我在测试中还注意到一个细节:Koopman MPC在靠近竖直向上平衡点时,稳态精度反而比工作点线性化MPC差一点点。原因是Koopman模型本身就是一个近似模型,在极小误差范围内,它的增益可能存在微小的偏移。如果系统对稳态精度要求极高,可以在平衡点附近切换到一个精确线性化模型,或者直接加一个积分环节来消除稳态误差。
5.2 字典维度、预测时域和采样周期的影响
参数影响这一块,我分别做了几组消融测试,结论比较有代表性:
- 字典维度:字典从8维加到15维时,预测误差下降非常明显;但从15维加到25维时,提升有限,甚至出现过拟合导致的验证集误差上升。对倒立摆这种系统,10到15个字典函数基本够用,没必要盲目加高维。
- 预测时域Np:Np太短(小于10)时,MPC看不到系统远期的动力学,控制效果明显变差;Np在20到40之间效果稳定;Np超过50之后,QP求解时间变长,但控制性能几乎不再提升。
- 采样周期dt:Koopman模型是离散模型,dt取得太大时欧拉近似的截断误差会比较大,K矩阵辨识出来也会失真。dt从0.02秒改成0.05秒后,控制精度明显下降。这说明Koopman MPC虽然让优化变简单了,并没有免除模型离散化的精度要求。
总的来说,调参顺序建议是:先固定采样周期,然后增加字典维度直到预测误差不再明显下降,最后再调Np和权重矩阵。不要一上来就同时调所有参数,那样很难定位问题。
5.3 从预测器到状态估计:Koopman线性空间的红利
标题里提到了“状态估计”,我在做这个Demo之前一直觉得Koopman只是换个方式做辨识,直到后来仔细想才意识到,线性预测器给状态估计带来的便利是实实在在的。
Koopman MPC在运行时,第一步就是把当前测量x0提升为z0=ψ(x0),然后用线性模型Az+Bu做未来预测。这个预测器本身就是一个线性递推结构。如果你需要做输出反馈控制,在Koopman提升空间里,你可以直接套用线性卡尔曼滤波或龙伯格观测器来估计提升状态,而不是在原始非线性系统上使用扩展卡尔曼滤波器去反复计算雅可比矩阵。
实际效果上,扩展卡尔曼滤波的雅可比计算在每一个采样周期都要做,而Koopman线性模型配合线性卡尔曼滤波只需要常数矩阵运算,计算量小得多,而且在字典函数设计得当的情况下,对强非线性系统的状态估计精度往往优于局部线性化的EKF。这一点对于想在嵌入式平台上做实时状态估计的人来说非常实用。
6. 我在实操中踩过的坑,以及我对Koopman MPC的判断
6.1 数据与归一化:最容易被忽视的一环
第一个坑是训练数据激励不充分。我有一次图省事,只在平衡点附近用小幅值随机扰动采集数据,结果Koopman模型在小偏差范围内表现很好,但MPC一遇到大偏差工况就失灵。后来我才意识到,EDMD本质上是在拟合一个全局模型,训练数据必须覆盖你想控制的整个状态空间范围。如果训练数据里根本没有出现过大角度摆动的数据点,Koopman矩阵自然学不会怎么处理大摆角。
第二个坑是量纲差异。倒立摆的状态里,角度(大约0到3)和角速度(可能达到每秒几弧度)数值量级差异很大;如果直接用原始数值做字典函数,回归时小量纲的状态可能被忽略。我的做法是先把状态归一化到相近的量级,再构造字典,拟合完A、B之后在MPC层面对提升状态本身进行反归一化处理。这一步虽然麻烦,但能显著提高数值稳定性。
第三个坑是参考点的常数项。如果你在字典里加了常数项1,那么z_ref的常数项也必须是1。我在第一次测试时把z_ref直接设成lift_fun([pi;0]),发现模型输出怎么都稳定不到目标点,查了很久才意识到是常数项偏置导致的问题。这个问题看起来小,实际影响巨大。
6.2 提升空间约束与可行性保护
提升空间是n_lift维的,但只有前nx个维度能直接对应原始状态。经常有人把约束施加在全部提升状态上,比如对sin(x₁)这个维度也加一个上下界,这完全没有物理意义,只会让QP可行域变得别扭甚至无解。
正确的做法是只对C·z里的原始状态部分施加约束。同时还要注意,预测模型毕竟是一个近似模型,预测时域拉得越长,模型偏差累积就越大。对于大范围工况,建议在MPC里加入一个“可行域保护”:如果求解出来的第一步控制量明显异常,比如超出输入限幅或者让状态预测值发散,控制器应该切换到更保守的控制策略,而不是盲目相信QP的结果。
我在这方面的处理是,实时监测模型预测的末端状态范数,一旦发现超出经验范围,就缩短预测时域重新求解。这个方法不算严谨,但在工程里很管用。
6.3 什么场景适合,什么场景别硬上
基于这个项目和一些模拟实验,我对Koopman MPC的适用场景形成了几点判断。
适合的场景包括:状态维度适中(2到10维左右),非线性比较光滑(没有硬切换、碰撞、死区这类不连续动态),而且工作范围比较明确,可以提前采集到覆盖整个范围的训练数据。机器人关节控制、连续搅拌反应釜的温度控制、飞行器姿态控制等都属于典型受益场景。
不适合的场景包括:系统含强不连续动态(比如齿轮间隙、摩擦切换)、状态维度很高但训练数据量严重不足、或者系统本身会频繁运行到训练覆盖范围之外的区域。在这些情况下,Koopman模型的近似误差会迅速放大,MPC的可靠性很难保证。
如果非要在这类场景里用,也有一条变通路线:引入鲁棒修正项,把Koopman模型的近似误差当作有界扰动,在MPC约束上做收缩;或者配合在线增量学习方法,每隔一段时间用最新数据重新拟合A、B矩阵。这些属于更进阶的方向,我自己还在探索中。
6.4 后续可以朝哪些方向扩展
Koopman MPC这套框架的扩展空间其实比想象中更大。我现在在尝试的方向主要有三个:
一是输出反馈合并。既然提升空间里有现成的线性预测器,下一步可以在这个线性状态空间模型上直接设计线性观测器或者卡尔曼滤波,解决系统只有部分状态可测的问题。二是鲁棒Koopman MPC。把EDMD辨识得到的残差建模为有界不确定性,设计tube MPC或者min-max MPC,用来抵御模型误差的影响。三是与强化学习结合。利用Koopman线性模型作为一个可微分的环境模型,帮助策略优化更快收敛,这一点在机器人运动控制里非常有吸引力。
最后再分享一个我的个人体会:如果你的目标是快速验证一个Koopman MPC算法,请先把数据质量、字典设计和模型验证这三件事做扎实,再去优化求解器参数。我踩过的最大的坑,不是quadprog调得不够好,而是Koopman模型本身在预测时域内就已经失真了,控制器再聪明也没用。先让预测器足够好,再谈控制器性能,这句话值得写在工位上。