MPC这东西圈内聊得最多的不是算法本身多高深而是“同一个MPC换一种模型写法效果能差出十万八千里”。尤其是输入增量形式工程里几乎绕不开——执行器有物理限制你总不希望控制量在每个采样周期里都剧烈跳变。但很多人刚开始接触时容易被各种“状态空间MPC公式”绕晕明明都是预测控制怎么一会用原始状态、一会用增广状态预测矩阵的形式还各不相同这篇就把这块彻底捋清楚我会用Matlab代码逐一实现不同状态空间公式对比它们的推导过程、实现细节和实际效果。适合刚入门MPC、想从代码层面吃透公式的研究生以及做运动控制、过程控制、机器人轨迹跟踪的工程师。1. MPC公式化的两条路线为什么“输入增量”是分水岭先从最朴素的问题开始MPC要做的事本质上是“根据当前状态预测未来一段时间的输出然后求解一个带约束的优化问题得到最优控制序列”。但“控制量”怎么进优化问题直接决定了公式长什么样。最常见的分水岭就是你优化的是绝对控制量 u(k), u(k1), ...还是优化的是控制增量 Δu(k), Δu(k1), ...。1.1 绝对控制量形式看着简单实际用起来处处受限最原始的状态空间MPC用的是x(k1) A x(k) B u(k)y(k) C x(k)这里优化变量是 u(k), ..., u(kN-1)N 是预测时域。目标函数里直接写 u 的二次项约束也直接约束 u 的上下限。听起来非常直接但问题在于工程上你几乎不可能只约束 u 而不关心它的变化率。电机有电流环限幅阀门有开启速率限制机械臂的关节力矩也不能瞬间从正最大跳到负最大。如果你只优化 u那求解器完全可能给出相邻两步之间剧烈抖动的控制序列而且这种抖动在目标函数里是“合法”的。当然你说我可以在目标函数里加一项 ||u(k) - u(k-1)||² 的惩罚这确实能抑制抖动但注意这一项用的是“上一步的 u(k-1)”而 u(k-1) 本身是过去的已知量它在当前优化问题里是常数。这种“半吊子”的增量惩罚会让问题变得不对称——你并没有把增量真正变成决策变量只是给变化量加了软惩罚约束层面依旧管不住速率。1.2 输入增量形式为什么是工业主流真正干净的做法是把决策变量直接换成 Δu(k), ..., Δu(kN-1)并且把模型改写成以“当前状态 上一时刻控制量”为扩展状态的形式。这样做的直接收益有三个第一约束能力覆盖全面。你既可以约束 u 的幅值也可以约束 Δu 的速率甚至可以同时约束输出 y 的上下限因为所有约束都可以写成 Δu 序列的线性组合。第二积分作用天然内嵌。增量形式的MPC本质上描述了“控制量本身也是被积分出来的状态”模型中的输入矩阵发生了改变闭环系统对恒定扰动有天然的抑制能力不需要额外加积分器。第三目标函数里的权重矩阵物理意义清晰。Q 对应输出误差R 对应“控制量变化的代价”而不是“控制量本身的代价”。在工程调参时你只需要问自己“我希望输出误差和控制变化速度之间怎么取舍”这比直接调“控制量绝对值”的权重要直观得多。所以我们要研究的几种不同状态空间MPC公式本质上是围绕“如何把增量形式写进状态空间模型”展开的。有人用非增广状态 重写预测方程有人用增广状态 直接套标准MPC还有人在此基础上再扩展扰动模型。每种的推导复杂度、数值特性、代码实现量都不同。2. 三种主流状态空间MPC公式的推导与对比这里我会把三种公式从头推导一遍。为了表述统一假设被控对象是线性离散时间系统x(k1) A x(k) B u(k)y(k) C x(k)其中 x ∈ R^nu ∈ R^my ∈ R^p。采样周期已经体现在 A、B 中了后面不再重复写时间标度。2.1 公式一非增广模型直接推导增量预测方程这种思路不改变状态的定义仍然用原始 x 作为状态但在构造预测方程时用 Δu 而不是 u 来表示未来的 x。做法是把 u(ki) 拆成u(ki) u(k-1) Σ_{j0}^{i} Δu(kj)代入状态递推式x(k1) A x(k) B [u(k-1) Δu(k)]x(k2) A x(k1) B [u(k-1) Δu(k) Δu(k1)] A² x(k) A B u(k-1) A B Δu(k) B Δu(k1) B Δu(k)以此类推N 步之后可以写成紧凑形式X Φ_x x(k) Γ_u u(k-1) Ψ ΔU其中 X [x(k1); x(k2); ...; x(kN)]ΔU [Δu(k); Δu(k1); ...; Δu(kN-1)]矩阵 Φ_x、Γ_u、Ψ 都是由 A、B 组合而成。输出预测是 Y C_x X其中 C_x 是分块对角矩阵每个对角块是 C。这个公式的好处是状态维度没有扩大保持 n 维计算预测矩阵时用的是原始 A、B看起来和标准MPC最接近。缺点也很明显你需要在预测方程里额外拎出来一个 u(k-1) 项推导时容易漏掉而且当你想把“输出约束”“控制幅值约束”统一写成矩阵形式时约束矩阵需要单独运算代码不如增广形式规整。2.2 公式二增广状态空间模型经典输入增量形式这是目前教材和工控软件里最流行的形式。定义增广状态ξ(k) [ x(k); u(k-1) ]那么ξ(k1) [ A B; 0_m×n I_m ] ξ(k) [ B; I_m ] Δu(k)记 A_tilde [A B; 0 I]B_tilde [B; I]输出方程y(k) [ C 0 ] ξ(k) C_tilde ξ(k)这样就把“增量”彻底融进了模型里从形式上看它和标准MPC公式完全一样只是 A、B、C 换成了增广版本。后续的预测矩阵、QP问题构建可以直接套用所有成熟的理论和代码模板。增广形式最大的优点就是“规整”。预测方程可以写成Y F ξ(k) G ΔU其中 F 和 G 的结构非常统一无论预测时域取多少只要用同样的递推规则就能生成。而且控制幅值约束写起来非常舒服因为 u(ki) u(k-1) Σ_{j0}^{i} Δu(kj)本身就是增广状态的一部分约束矩阵可以写成分块三角矩阵Matlab里用几个 tril、ones 就能拼出来。代价是状态维度变成 nm。如果原始系统有4个状态、2个输入增广后就变成6维计算量会略微增加但对现代计算机和嵌入式处理器来说只要系统规模不是几百上千维这点代价完全可以忽略。2.3 公式三带扰动模型的增广MPC无偏跟踪的关键第三种公式在高性能控制里经常被提及就是在增广状态里再加入“扰动状态”。模型的来源通常是觉得“模型和真实系统之间总存在失配被控过程还会有外来扰动”。为了在MPC框架内对抗这种恒值扰动构造如下增广模型x(k1) A x(k) B u(k) B_d d(k)d(k1) d(k)把扰动 d 当成一个状态它的动态是恒值。然后定义新的增广状态ξ_d(k) [ x(k); u(k-1); d(k) ]这样一步到位得到更高维的模型。输出预测中扰动项会直接体现在预测值里MPC优化时会主动产生一个“抵消扰动”的控制量从而实现无偏跟踪。当然这要求扰动可观测不然你只知道它的存在却不知道它多大。在实际代码中通常会用状态观测器比如卡尔曼滤波去估计 d或者在简单情况下用输出误差的动态来隐式估计。这个公式的适用场景是系统存在模型失配、恒定负载扰动且你希望稳态误差为零。比如电机控制里的负载转矩突变、温度控制里的环境温度变化。缺点是状态维度进一步增加观测器需要单独设计代码复杂度显著上升。2.4 三种公式的快速对比公式状态定义决策变量是否带积分作用代码复杂度适用场景非增广直接推导原始 xΔu有通过增量项中等预测矩阵推导繁琐教学演示、系统维度极低时经典增广[x; u(k-1)]Δu有低结构规整绝大多数线性MPC工程应用带扰动增广[x; u(k-1); d]Δu强扰动补偿高需观测器高精度伺服、抗扰控制这个表只是概览实际选型需要结合模型复杂度、约束数量和执行周期。下面我重点讲经典增广的Matlab实现因为它是三个公式里性价比最高的一个。3. Matlab代码实现经典增广MPC的完整拆解3.1 问题设定与离散化假设我们研究一个典型的二阶系统比如直流电机的简化模型状态是电流和转速控制量是电压输出是转速。连续时间状态空间矩阵给出如下这里用虚构参数方便复现A_c [0 -10; 1 -2]B_c [10; 0]C_c [0 1]D_c 0取采样周期 Ts 0.05s用 Matlab 的 c2d 做零阶保持离散化A_c [0 -10; 1 -2]; B_c [10; 0]; C_c [0 1]; D_c 0; Ts 0.05; sys_c ss(A_c, B_c, C_c, D_c); sys_d c2d(sys_c, Ts, zoh); A sys_d.A; B sys_d.B; C sys_d.C;这里的 C 是行向量输出转速 y 是标量所以 p 1n 2m 1。系统很简单但足够说明问题。3.2 构建增广状态空间矩阵按第二节的公式二增广状态是 [x; u(k-1)][n, ~] size(A); [m, ~] size(B); A_tilde [A, B; zeros(m, n), eye(m)]; B_tilde [B; eye(m)]; C_tilde [C, zeros(size(C,1), m)];这一步千万别写错。很多新手容易把右下角写成零矩阵那 u(k-1) 就会在状态更新时被抹掉系统变成“伪积分”结构闭环可能不稳定。写成 eye(m) 表示 u(k-1) 会原封不动保留到下一时刻直到被新的 Δu 覆盖。3.3 预测矩阵 F 和 G 的生成预测时域设为 N 10控制时域 M 3通常 M N表示后续控制增量假设为零。我们需要两个矩阵Y_pred F * xi(k) G * dU其中 xi(k) 是当前增广状态dU 是 M 个未来控制增量。F 的大小是 (pN) × (nm)G 的大小是 (pN) × (m*M)。生成方法有两种循环累乘和矩阵幂。我推荐循环不容易出错N 10; M 3; F zeros(N, nm); G zeros(N, M); % 本质上是模拟从当前状态开始逐步递推 % F 的第 i 块对应 C_tilde * A_tilde^(i1) % G 的第 i 行是 C_tilde * [A_tilde^(i) * B_tilde, A_tilde^(i-1)*B_tilde, ...] for i 1:N F(i, :) C_tilde * A_tilde^i; for j 1:M if j i G(i, j) C_tilde * A_tilde^(i-j) * B_tilde; end end end这里 G(i,j) 的含义是第 i 步预测输出对第 j 个控制增量的敏感系数。注意索引从1开始而数学递推中 k1 是第一步。如果预测时域从 k1 开始那么 G 的第一行应该是 C_tilde * B_tilde对应 Δu(k) 对 y(k1) 的直接影响。上面的代码用 A_tilde^0 I 来算当 ij1 时A_tilde^(1-1) * B_tilde B_tilde正好对。3.4 目标函数与约束的矩阵化经典增广形式的目标函数写成J (Y_ref - Y_pred) Q (Y_ref - Y_pred) dU R dU其中 Q 是输出误差权重对角分块R 是控制增量权重对角分块。为了转成标准二次规划QP形式min 0.5 * dU H dU f dU我们展开 J丢掉常数项Y_ref Q Y_ref 和交叉项中不包含 dU 的常数部分得到H 2 * (G Q G R)f -2 * G Q (Y_ref - F * xi_k)注意这里用 2 是因为 quadprog 求解的是 1/2 dU H dU f dU所以 H 要乘以 2。如果不乘结果也不会错太大但目标函数与约束的配比就不是你设计的那样了。Q 10 * eye(N); R 0.5 * eye(M); H 2 * (G * Q * G R);约束包含三块第一控制增量约束lb_dU dU ub_dU。直接作为优化变量的上下界传入 quadprog。第二控制幅值约束:u_min u(k-1) 累积增量 u_max。写起来是线性不等式约束 A_dU * dU b_dU其中累积增量矩阵是下三角矩阵。第三输出约束y_min Y_pred y_max。由于 Y_pred F xi_k G dU所以也是线性不等式只是需要把 F xi_k 这一项移到右侧。u_min -1; u_max 1; du_min -0.2; du_max 0.2; % 控制幅值约束 % u(ki) u(k-1) sum_{j0}^{i} du(kj) % 写成矩阵形式tril(ones(M,M)) * du u_k_1 u_min and u_max Tril tril(ones(M, M)); % 下三角累积矩阵 % 注意这里的预测步数是 M因为只有 M 个增量变量可调后续增量为零 % 但我们的预测矩阵 G 用的是 N 步预测控制量在 M 步以后保持不变 A_cons [Tril; -Tril]; b_cons [u_max - u_k_1; -(u_min - u_k_1)];上面的 u_k_1 是上一时刻的实际控制量在每一拍都是已知标量。3.5 QP求解与滚动时域用 quadprog 求解每个采样周期的问题options optimoptions(quadprog, Display, off, Algorithm, active-set); [dU_opt, ~, exitflag] quadprog(H, f, A_cons, b_cons, [], [], du_min*ones(M,1), du_max*ones(M,1), [], options);注意 quadprog 要求 H 是半正定所以 R 矩阵必须大于等于零。实际调试时 R 取太小会导致 H 接近奇异求解器会报警或解出离谱的值这也是后面常见问题要重点提的。滚动时域的意思是只取 dU_opt(1) 作为当前时刻的控制增量更新 u(k) u(k-1) dU_opt(1)然后施加给系统到下一采样周期再重复整个求解。代码框架xi [x0; u_prev]; u u_prev; for k 1:sim_steps % 测量/获取当前状态 x % 计算 Y_ref参考轨迹可以预知或实时给 % 构造 F、G、H、f、约束 % 求解 dU_opt dU dU_opt(1); u u_prev dU; % 用 u 驱动真实系统仿真x_next A*x B*u % 更新 u_prev u % 更新 xi [x; u_prev] end这里有个关键点在仿真中下一步状态要用真实模型或者被控对象提供的测量值来算不要用增广模型里的 A_tilde因为增广模型里的 u(k-1) 只是状态拍真实物理量是 x 和 u。3.6 完整闭环仿真的范例我给出一个可以直接跑的完整脚本注释写得比较细% 参数 Ts 0.05; N 10; M 3; Q 10 * eye(N); R 0.5 * eye(M); u_min -1; u_max 1; du_min -0.2; du_max 0.2; % 模型虚构的二阶系统 A_c [0 -10; 1 -2]; B_c [10; 0]; C_c [0 1]; D_c 0; sys_c ss(A_c, B_c, C_c, D_c); sys_d c2d(sys_c, Ts, zoh); A sys_d.A; B sys_d.B; C sys_d.C; % 增广 [n, ~] size(A); [m, ~] size(B); A_tilde [A, B; zeros(m,n), eye(m)]; B_tilde [B; eye(m)]; C_tilde [C, zeros(1,m)]; % 预计算F,G对于时不变系统可以只算一次但为了直观这里在循环外算 F zeros(N, nm); G zeros(N, M); for i 1:N F(i, :) C_tilde * A_tilde^i; for j 1:M if j i G(i, j) C_tilde * A_tilde^(i-j) * B_tilde; end end end % 仿真设置 sim_steps 200; x [0; 0]; u_prev 0; y_ref 1; % 单位阶跃 Y_ref y_ref * ones(N, 1); history_x zeros(n, sim_steps); history_u zeros(1, sim_steps); history_y zeros(1, sim_steps); for k 1:sim_steps xi [x; u_prev]; % 目标函数 H 2 * (G * Q * G R); f -2 * G * Q * (Y_ref - F * xi); % 约束 A_cons [tril(ones(M,M)); -tril(ones(M,M))]; b_cons [u_max - u_prev; -(u_min - u_prev)]; % 求解 options optimoptions(quadprog, Display, off, Algorithm, active-set); [dU_opt, ~, ~] quadprog(H, f, A_cons, b_cons, [], [], ... du_min*ones(M,1), du_max*ones(M,1), [], options); dU dU_opt(1); u u_prev dU; % 系统动态 x A * x B * u; y C * x; u_prev u; % 记录 history_x(:, k) x; history_u(k) u; history_y(k) y; end % 绘图 t (0:sim_steps-1) * Ts; figure; subplot(2,1,1); plot(t, history_y, LineWidth, 1.5); hold on; plot(t, ones(1,sim_steps), --); xlabel(时间/s); ylabel(输出); subplot(2,1,2); stairs(t, history_u, LineWidth, 1.5); xlabel(时间/s); ylabel(控制量);这段代码跑下来你应该能看到输出大约经过几次调整后稳定到 1控制量被限制在 [-1,1] 之间而且变化速率被限在 ±0.2。如果没有速率约束控制量可能会一步跳满这就是增量约束的价值。4. 关键体会三种公式切换时的坑与Matlab调试技巧4.1 预测矩阵维度不匹配是最常见的问题我见过太多人把增广状态搞对之后预测矩阵 F 和 G 的维度跟 C_tilde 相乘出错。你只要记得一句话F 的每一行对应一个预测步行的长度等于增广状态的维度G 的每一行也对应一个预测步列数等于控制增量个数通常就是 M。如果 F 是 N×(nm)G 是 N×M那 Y_pred 就是 N×1Y_ref 也必须是 N×1。检查维度时直接用 size 打印出来逐一对比盯着公式猜要快得多。4.2 控制时域 M 小于预测时域 N 的矩阵生成秘技如果 M N请注意 G 矩阵里只有前 M 列有值从第 M1 列开始应该全为零。上面的代码里内层 j 循环只到 M且 if j i 保证了当 i 较小时j 太大对应的项不存在。这里有个容易忽略的点当 i 很大、j 很小的时候A_tilde^(i-j) 的幂次很大数值上可能会出现近乎奇异的项这时应该考虑用递推而不是矩阵幂生成 G。对于一般的稳定系统N 不超过 20 问题不大如果 N 很大建议把生成 G 的循环改写成递推累加避免大幂次导致精度损失。4.3 quadprog 的算法选择与目标函数缩放Matlab 的 quadprog 默认算法在不同版本下可能不同有些老版本默认 active-set新版本默认 interior-point-convex。对于 MPC 这个规模不大、约束线性且上下界简单的 QP两者都能解。但注意 interior-point 算法可能会在约束比较紧时给出边界附近的解而 active-set 对初始点敏感。我通常的做法是显式指定 active-set因为 MPC 的 QP 一般只有几百个决策变量active-set 迭代次数少行为也容易预测。另外如果 H 的条件数太大因为 Q 和 R 相差几个数量级建议先对 Q 和 R 做归一化比如都除以输出量纲的平方否则求解器容易报“H 必须是正定”的警告。4.4 参考轨迹的预知与实时性问题上面代码里 Y_ref 全程是常数这只是最简化。如果参考轨迹是时变的需要在每个采样周期更新 Y_ref 向量。对于可以预知的参考轨迹比如规划的轨迹可以直接把未来 N 步的参考值放进去这就是“前馈式”MPC能显著改善跟踪滞后。如果你的应用只能拿到当前参考值那就让 Y_ref y_ref(k) * ones(N,1)也能用只是跟踪会有纯滞后。这在运动控制里要特别注意因为轨迹跟踪的滞后往往就是追踪误差的主要来源。4.5 三种公式在Matlab里的实际调试体验非增广公式在实现时最大的坑在于构造 Ψ 矩阵时容易把 u(k-1) 的传播项搞错。我曾经在调试时漏了 Γ_u 那一项结果系统响应奇慢后来发现是少加了一项 B*u(k-1) 的累加项。增广公式就完全不会出这种问题因为状态里带着 u(k-1)所有项都在预测矩阵里自动处理了。带扰动的公式调试难度主要在观测器如果 d 的估计滞后MPC 的补偿作用会来得慢这时需要检查扰动模型的收敛速度。4.6 给新手的Matlab调试建议我自己的习惯是三步走第一步先把未加约束的MPC闭环仿真跑通。只保留目标函数不加入任何不等式约束。如果这一步都不稳定那就是模型或预测矩阵错了此时去检查 A_tilde 和 B_tilde。第二步加入控制幅值约束。观察控制量是否触碰边界触碰边界时闭环是否仍稳定。如果触碰边界后发散多半是约束矩阵 A_cons 写错了尤其是符号问题。第三步加入速率约束和控制时域缩减。速率约束是增量形式的专属优势加上后观察响应变慢的程度再调 R。每加一层约束都单独验证千万不要一上来全加上否则出了问题根本分不清是预测错误还是约束冲突。5. 参数整定的经验笔记与常见问题速查5.1 预测时域N、控制时域M和采样周期Ts怎么配先说采样周期它必须是系统主要动力学的 1/5 ~ 1/10。比如系统开环时间常数是 1s那采样周期取 0.1 ~ 0.2s 比较稳妥。Ts 太大会丢失两采样点之间的动态太小会让每一步控制增量都很小导致跟踪迟钝。预测时域 N 的物理含义是“看多远的未来”至少要覆盖系统动态的主要时间跨度。一个简单经验N * Ts ≈ 3 ~ 5 倍系统时间常数。比如时间常数 1sTs0.1s那 N 至少取 30 到 50。取小了MPC 只看近处约束处理和稳定性都会变差。控制时域 M 决定了“自由度的数量”。M 越大控制越激进计算量也越大M 很小比如1控制量几乎只在第一拍决定后面全不更新响应会很保守。一般取 N 的 1/3 到 1/2 作为 M。5.2 Q和R的整定顺序跟 PID 调参一样先调 Q 让跟踪误差尽快收敛再用 R 抑制控制量的高频抖振。但要注意增量形式的 R 直接惩罚的是 Δu 的幅值R 太小控制增量会很大但因为有速率约束du_min/du_max兜底一般不会失控只是会让 u 频繁触及约束边界。R 太大系统响应变慢而且容易产生一个现象叫“提前减速”——MPC 预测到未来误差不大所以提前开始减小控制增量导致跟踪末端有明显的拖尾。遇到这种情况把 R 减半试试。5.3 约束冲突怎么办当输出约束、控制幅值约束和控制速率约束同时存在时可能出现“无可行解”。最简单的手段是扩大控制时域 M或者放宽速率约束。另一种更工程化的做法是将输出约束转成软约束引入松弛变量在目标函数里加一个惩罚项。Matlab里实现软约束比较简单你只需要把原问题里 A_cons 对应的输出约束改成 A_cons * dU b_cons slack同时在目标函数里加上 lambda * slack。这里 slack 是额外的决策变量H 矩阵对应位置为零f 向量对应位置加 lambda。这样无论约束是否冲突求解器都能给出一个解只是当约束被突破时slack 变大目标函数上升表现为“尽量满足约束但允许在极端情况下微量突破”。5.4 常见问题速查表现象可能原因排查方法闭环发散A_tilde 或 B_tilde 构造错误打印闭环极点检查增广矩阵的特征值仿真速度极慢每次循环里重复构造 F、G时不变系统 F、G 移到循环外输出稳态误差不为零模型无积分作用或没有扰动补偿改用带扰动增广公式或检查权重 Q控制量剧烈振荡R 太小或速率约束太松增大 R或收紧 du_min/du_maxquadprog 报错“H 不正定”R 为零或 Q 与 R 数量级悬殊给 R 加一个很小的正则项比如 1e-6参考轨迹跟踪滞后明显Y_ref 没有使用未来值用预知的参考轨迹填充未来 N 步控制量始终卡在下界约束写反了符号检查 A_cons 和 b_cons 的正负号5.5 一个调参案例我实际跑过的一个例子模拟项目X是一个温度控制系统采样周期0.5s系统时间常数约10s。最初取 N20M5Q1R0.1响应非常慢因为 N*Ts10s只覆盖了1倍时间常数MPC 看不到更远的未来。后来把 N 加到 40覆盖20sM 加到 10R 降到 0.02输出响应明显加快最终稳态误差为零。这就是“预测时域必须覆盖系统主要动态”的直观教训。还有一次我把速率约束设成 du_max0.5而实际系统每0.5s采样一次相当于控制量每秒变化不超过1这本身是合理的。但 R 设成了 0.001结果每步控制增量都碰到约束边界控制量一路滑到底系统超调接近30%。后来把 R 提到 0.05超调降到 8%。这说明速率约束和权重 R 相辅相成约束决定硬边界R 决定在边界内部的软取舍。6. 带扰动公式的Matlab扩展要点前面提到了第三种公式带扰动增广这里再补充一点代码层面的细节。扰动 d 通常不是直接可测的你需要用观测器估计。最简单的是一维扰动取扩展状态x_aug [x; u_prev; d]观测器设计可以用标准的扩展卡尔曼滤波器也可以用一个简单的龙伯格观测器前提是系统可观。Matlab里可以这样构造观测矩阵% 假设扰动模型: d(k1) d(k) % 观测方程对输出 y 的扰动输入是 C_d通常等于 [C, 0, 0; 0, 0, 1] 的形式 % 具体矩阵根据你的对象模型来定在MPC预测矩阵里F 和 G 需要对应扩展后的状态维度所以 F 的列为 n1_dim_u1_dim_dG 的列还是 M。扰动状态在预测中起到的作用是预测输出 名义模型输出 扰动在输出上的投影这样预测误差不再只是“当前误差”而是被扰动模型解释并持续到未来。实际调试中我建议先把扰动公式用在无扰动环境里验证确认闭环性能与经典增广完全一致再人为加入阶跃扰动观察MPC是否有补偿能力。如果发现补偿速度慢可以增大观测器增益但注意观测器增益过大会放大噪声这是一个需要权衡的点。7. 结束前的几个实操心得最后说点平时文档里不爱写的经验。第一做MPC仿真时千万别只盯着输出曲线和控制量曲线一定要把控制增量曲线也画出来。很多问题单看控制量觉得挺正常但控制增量已经高频抖动得不成样子这往往是 R 权重和约束不匹配的信号。我在做某图像处理Demo的辅助控制的时候就是靠画 Δu 曲线才发现一个采样周期里的数值振荡问题。第二代码模块化很重要。我说的是把“模型离散化”、“增广矩阵构造”、“预测矩阵生成”、“QP求解”、“仿真循环”分别写成独立的函数。初期调试用脚本没问题但一旦要换模型、换参数函数化能帮你节省大量时间。我自己的习惯是把预测矩阵生成写成 generatePredictionMatrices(A_tilde, B_tilde, C_tilde, N, M)输入输出一目了然后续切换公式时只需要替换这个函数内部实现。第三如果目标平台是嵌入式C代码Matlab里仿真用的 quadprog 通常不能直接移植。但你可以先用Matlab验证控制律然后把 QP 求解换成嵌入式实时求解器比如基于内点法或活动集法的轻量库。这个过程中增广公式的优势就体现出来了预测矩阵和约束矩阵的形式简单移植时不容易出错。相比之下非增广公式在移植时还得单独处理 u(k-1) 项容易埋bug。谁用谁知道。MPC没有想象中那么神秘但也没有某些论文写得那么轻巧。把输入增量这个点吃透不同状态空间公式对你来说就只是“怎么组织模型”的选择题而不是拦路虎。希望这篇能让你少走点弯路动手把代码跑起来之后你会发现MPC真正的乐趣在于“约束下求最优”的那个瞬间——一切看起来都在边界上跳舞但舞步是稳的。 
长河君 · 2026/10/14 6:55:48 · 阅读 1.2k MPC状态空间公式详解：输入增量形式与Matlab实现对比
MPC这东西圈内聊得最多的不是算法本身多高深而是“同一个MPC换一种模型写法效果能差出十万八千里”。尤其是输入增量形式工程里几乎绕不开——执行器有物理限制你总不希望控制量在每个采样周期里都剧烈跳变。但很多人刚开始接触时容易被各种“状态空间MPC公式”绕晕明明都是预测控制怎么一会用原始状态、一会用增广状态预测矩阵的形式还各不相同这篇就把这块彻底捋清楚我会用Matlab代码逐一实现不同状态空间公式对比它们的推导过程、实现细节和实际效果。适合刚入门MPC、想从代码层面吃透公式的研究生以及做运动控制、过程控制、机器人轨迹跟踪的工程师。1. MPC公式化的两条路线为什么“输入增量”是分水岭先从最朴素的问题开始MPC要做的事本质上是“根据当前状态预测未来一段时间的输出然后求解一个带约束的优化问题得到最优控制序列”。但“控制量”怎么进优化问题直接决定了公式长什么样。最常见的分水岭就是你优化的是绝对控制量 u(k), u(k1), ...还是优化的是控制增量 Δu(k), Δu(k1), ...。1.1 绝对控制量形式看着简单实际用起来处处受限最原始的状态空间MPC用的是x(k1) A x(k) B u(k)y(k) C x(k)这里优化变量是 u(k), ..., u(kN-1)N 是预测时域。目标函数里直接写 u 的二次项约束也直接约束 u 的上下限。听起来非常直接但问题在于工程上你几乎不可能只约束 u 而不关心它的变化率。电机有电流环限幅阀门有开启速率限制机械臂的关节力矩也不能瞬间从正最大跳到负最大。如果你只优化 u那求解器完全可能给出相邻两步之间剧烈抖动的控制序列而且这种抖动在目标函数里是“合法”的。当然你说我可以在目标函数里加一项 ||u(k) - u(k-1)||² 的惩罚这确实能抑制抖动但注意这一项用的是“上一步的 u(k-1)”而 u(k-1) 本身是过去的已知量它在当前优化问题里是常数。这种“半吊子”的增量惩罚会让问题变得不对称——你并没有把增量真正变成决策变量只是给变化量加了软惩罚约束层面依旧管不住速率。1.2 输入增量形式为什么是工业主流真正干净的做法是把决策变量直接换成 Δu(k), ..., Δu(kN-1)并且把模型改写成以“当前状态 上一时刻控制量”为扩展状态的形式。这样做的直接收益有三个第一约束能力覆盖全面。你既可以约束 u 的幅值也可以约束 Δu 的速率甚至可以同时约束输出 y 的上下限因为所有约束都可以写成 Δu 序列的线性组合。第二积分作用天然内嵌。增量形式的MPC本质上描述了“控制量本身也是被积分出来的状态”模型中的输入矩阵发生了改变闭环系统对恒定扰动有天然的抑制能力不需要额外加积分器。第三目标函数里的权重矩阵物理意义清晰。Q 对应输出误差R 对应“控制量变化的代价”而不是“控制量本身的代价”。在工程调参时你只需要问自己“我希望输出误差和控制变化速度之间怎么取舍”这比直接调“控制量绝对值”的权重要直观得多。所以我们要研究的几种不同状态空间MPC公式本质上是围绕“如何把增量形式写进状态空间模型”展开的。有人用非增广状态 重写预测方程有人用增广状态 直接套标准MPC还有人在此基础上再扩展扰动模型。每种的推导复杂度、数值特性、代码实现量都不同。2. 三种主流状态空间MPC公式的推导与对比这里我会把三种公式从头推导一遍。为了表述统一假设被控对象是线性离散时间系统x(k1) A x(k) B u(k)y(k) C x(k)其中 x ∈ R^nu ∈ R^my ∈ R^p。采样周期已经体现在 A、B 中了后面不再重复写时间标度。2.1 公式一非增广模型直接推导增量预测方程这种思路不改变状态的定义仍然用原始 x 作为状态但在构造预测方程时用 Δu 而不是 u 来表示未来的 x。做法是把 u(ki) 拆成u(ki) u(k-1) Σ_{j0}^{i} Δu(kj)代入状态递推式x(k1) A x(k) B [u(k-1) Δu(k)]x(k2) A x(k1) B [u(k-1) Δu(k) Δu(k1)] A² x(k) A B u(k-1) A B Δu(k) B Δu(k1) B Δu(k)以此类推N 步之后可以写成紧凑形式X Φ_x x(k) Γ_u u(k-1) Ψ ΔU其中 X [x(k1); x(k2); ...; x(kN)]ΔU [Δu(k); Δu(k1); ...; Δu(kN-1)]矩阵 Φ_x、Γ_u、Ψ 都是由 A、B 组合而成。输出预测是 Y C_x X其中 C_x 是分块对角矩阵每个对角块是 C。这个公式的好处是状态维度没有扩大保持 n 维计算预测矩阵时用的是原始 A、B看起来和标准MPC最接近。缺点也很明显你需要在预测方程里额外拎出来一个 u(k-1) 项推导时容易漏掉而且当你想把“输出约束”“控制幅值约束”统一写成矩阵形式时约束矩阵需要单独运算代码不如增广形式规整。2.2 公式二增广状态空间模型经典输入增量形式这是目前教材和工控软件里最流行的形式。定义增广状态ξ(k) [ x(k); u(k-1) ]那么ξ(k1) [ A B; 0_m×n I_m ] ξ(k) [ B; I_m ] Δu(k)记 A_tilde [A B; 0 I]B_tilde [B; I]输出方程y(k) [ C 0 ] ξ(k) C_tilde ξ(k)这样就把“增量”彻底融进了模型里从形式上看它和标准MPC公式完全一样只是 A、B、C 换成了增广版本。后续的预测矩阵、QP问题构建可以直接套用所有成熟的理论和代码模板。增广形式最大的优点就是“规整”。预测方程可以写成Y F ξ(k) G ΔU其中 F 和 G 的结构非常统一无论预测时域取多少只要用同样的递推规则就能生成。而且控制幅值约束写起来非常舒服因为 u(ki) u(k-1) Σ_{j0}^{i} Δu(kj)本身就是增广状态的一部分约束矩阵可以写成分块三角矩阵Matlab里用几个 tril、ones 就能拼出来。代价是状态维度变成 nm。如果原始系统有4个状态、2个输入增广后就变成6维计算量会略微增加但对现代计算机和嵌入式处理器来说只要系统规模不是几百上千维这点代价完全可以忽略。2.3 公式三带扰动模型的增广MPC无偏跟踪的关键第三种公式在高性能控制里经常被提及就是在增广状态里再加入“扰动状态”。模型的来源通常是觉得“模型和真实系统之间总存在失配被控过程还会有外来扰动”。为了在MPC框架内对抗这种恒值扰动构造如下增广模型x(k1) A x(k) B u(k) B_d d(k)d(k1) d(k)把扰动 d 当成一个状态它的动态是恒值。然后定义新的增广状态ξ_d(k) [ x(k); u(k-1); d(k) ]这样一步到位得到更高维的模型。输出预测中扰动项会直接体现在预测值里MPC优化时会主动产生一个“抵消扰动”的控制量从而实现无偏跟踪。当然这要求扰动可观测不然你只知道它的存在却不知道它多大。在实际代码中通常会用状态观测器比如卡尔曼滤波去估计 d或者在简单情况下用输出误差的动态来隐式估计。这个公式的适用场景是系统存在模型失配、恒定负载扰动且你希望稳态误差为零。比如电机控制里的负载转矩突变、温度控制里的环境温度变化。缺点是状态维度进一步增加观测器需要单独设计代码复杂度显著上升。2.4 三种公式的快速对比公式状态定义决策变量是否带积分作用代码复杂度适用场景非增广直接推导原始 xΔu有通过增量项中等预测矩阵推导繁琐教学演示、系统维度极低时经典增广[x; u(k-1)]Δu有低结构规整绝大多数线性MPC工程应用带扰动增广[x; u(k-1); d]Δu强扰动补偿高需观测器高精度伺服、抗扰控制这个表只是概览实际选型需要结合模型复杂度、约束数量和执行周期。下面我重点讲经典增广的Matlab实现因为它是三个公式里性价比最高的一个。3. Matlab代码实现经典增广MPC的完整拆解3.1 问题设定与离散化假设我们研究一个典型的二阶系统比如直流电机的简化模型状态是电流和转速控制量是电压输出是转速。连续时间状态空间矩阵给出如下这里用虚构参数方便复现A_c [0 -10; 1 -2]B_c [10; 0]C_c [0 1]D_c 0取采样周期 Ts 0.05s用 Matlab 的 c2d 做零阶保持离散化A_c [0 -10; 1 -2]; B_c [10; 0]; C_c [0 1]; D_c 0; Ts 0.05; sys_c ss(A_c, B_c, C_c, D_c); sys_d c2d(sys_c, Ts, zoh); A sys_d.A; B sys_d.B; C sys_d.C;这里的 C 是行向量输出转速 y 是标量所以 p 1n 2m 1。系统很简单但足够说明问题。3.2 构建增广状态空间矩阵按第二节的公式二增广状态是 [x; u(k-1)][n, ~] size(A); [m, ~] size(B); A_tilde [A, B; zeros(m, n), eye(m)]; B_tilde [B; eye(m)]; C_tilde [C, zeros(size(C,1), m)];这一步千万别写错。很多新手容易把右下角写成零矩阵那 u(k-1) 就会在状态更新时被抹掉系统变成“伪积分”结构闭环可能不稳定。写成 eye(m) 表示 u(k-1) 会原封不动保留到下一时刻直到被新的 Δu 覆盖。3.3 预测矩阵 F 和 G 的生成预测时域设为 N 10控制时域 M 3通常 M N表示后续控制增量假设为零。我们需要两个矩阵Y_pred F * xi(k) G * dU其中 xi(k) 是当前增广状态dU 是 M 个未来控制增量。F 的大小是 (pN) × (nm)G 的大小是 (pN) × (m*M)。生成方法有两种循环累乘和矩阵幂。我推荐循环不容易出错N 10; M 3; F zeros(N, nm); G zeros(N, M); % 本质上是模拟从当前状态开始逐步递推 % F 的第 i 块对应 C_tilde * A_tilde^(i1) % G 的第 i 行是 C_tilde * [A_tilde^(i) * B_tilde, A_tilde^(i-1)*B_tilde, ...] for i 1:N F(i, :) C_tilde * A_tilde^i; for j 1:M if j i G(i, j) C_tilde * A_tilde^(i-j) * B_tilde; end end end这里 G(i,j) 的含义是第 i 步预测输出对第 j 个控制增量的敏感系数。注意索引从1开始而数学递推中 k1 是第一步。如果预测时域从 k1 开始那么 G 的第一行应该是 C_tilde * B_tilde对应 Δu(k) 对 y(k1) 的直接影响。上面的代码用 A_tilde^0 I 来算当 ij1 时A_tilde^(1-1) * B_tilde B_tilde正好对。3.4 目标函数与约束的矩阵化经典增广形式的目标函数写成J (Y_ref - Y_pred) Q (Y_ref - Y_pred) dU R dU其中 Q 是输出误差权重对角分块R 是控制增量权重对角分块。为了转成标准二次规划QP形式min 0.5 * dU H dU f dU我们展开 J丢掉常数项Y_ref Q Y_ref 和交叉项中不包含 dU 的常数部分得到H 2 * (G Q G R)f -2 * G Q (Y_ref - F * xi_k)注意这里用 2 是因为 quadprog 求解的是 1/2 dU H dU f dU所以 H 要乘以 2。如果不乘结果也不会错太大但目标函数与约束的配比就不是你设计的那样了。Q 10 * eye(N); R 0.5 * eye(M); H 2 * (G * Q * G R);约束包含三块第一控制增量约束lb_dU dU ub_dU。直接作为优化变量的上下界传入 quadprog。第二控制幅值约束:u_min u(k-1) 累积增量 u_max。写起来是线性不等式约束 A_dU * dU b_dU其中累积增量矩阵是下三角矩阵。第三输出约束y_min Y_pred y_max。由于 Y_pred F xi_k G dU所以也是线性不等式只是需要把 F xi_k 这一项移到右侧。u_min -1; u_max 1; du_min -0.2; du_max 0.2; % 控制幅值约束 % u(ki) u(k-1) sum_{j0}^{i} du(kj) % 写成矩阵形式tril(ones(M,M)) * du u_k_1 u_min and u_max Tril tril(ones(M, M)); % 下三角累积矩阵 % 注意这里的预测步数是 M因为只有 M 个增量变量可调后续增量为零 % 但我们的预测矩阵 G 用的是 N 步预测控制量在 M 步以后保持不变 A_cons [Tril; -Tril]; b_cons [u_max - u_k_1; -(u_min - u_k_1)];上面的 u_k_1 是上一时刻的实际控制量在每一拍都是已知标量。3.5 QP求解与滚动时域用 quadprog 求解每个采样周期的问题options optimoptions(quadprog, Display, off, Algorithm, active-set); [dU_opt, ~, exitflag] quadprog(H, f, A_cons, b_cons, [], [], du_min*ones(M,1), du_max*ones(M,1), [], options);注意 quadprog 要求 H 是半正定所以 R 矩阵必须大于等于零。实际调试时 R 取太小会导致 H 接近奇异求解器会报警或解出离谱的值这也是后面常见问题要重点提的。滚动时域的意思是只取 dU_opt(1) 作为当前时刻的控制增量更新 u(k) u(k-1) dU_opt(1)然后施加给系统到下一采样周期再重复整个求解。代码框架xi [x0; u_prev]; u u_prev; for k 1:sim_steps % 测量/获取当前状态 x % 计算 Y_ref参考轨迹可以预知或实时给 % 构造 F、G、H、f、约束 % 求解 dU_opt dU dU_opt(1); u u_prev dU; % 用 u 驱动真实系统仿真x_next A*x B*u % 更新 u_prev u % 更新 xi [x; u_prev] end这里有个关键点在仿真中下一步状态要用真实模型或者被控对象提供的测量值来算不要用增广模型里的 A_tilde因为增广模型里的 u(k-1) 只是状态拍真实物理量是 x 和 u。3.6 完整闭环仿真的范例我给出一个可以直接跑的完整脚本注释写得比较细% 参数 Ts 0.05; N 10; M 3; Q 10 * eye(N); R 0.5 * eye(M); u_min -1; u_max 1; du_min -0.2; du_max 0.2; % 模型虚构的二阶系统 A_c [0 -10; 1 -2]; B_c [10; 0]; C_c [0 1]; D_c 0; sys_c ss(A_c, B_c, C_c, D_c); sys_d c2d(sys_c, Ts, zoh); A sys_d.A; B sys_d.B; C sys_d.C; % 增广 [n, ~] size(A); [m, ~] size(B); A_tilde [A, B; zeros(m,n), eye(m)]; B_tilde [B; eye(m)]; C_tilde [C, zeros(1,m)]; % 预计算F,G对于时不变系统可以只算一次但为了直观这里在循环外算 F zeros(N, nm); G zeros(N, M); for i 1:N F(i, :) C_tilde * A_tilde^i; for j 1:M if j i G(i, j) C_tilde * A_tilde^(i-j) * B_tilde; end end end % 仿真设置 sim_steps 200; x [0; 0]; u_prev 0; y_ref 1; % 单位阶跃 Y_ref y_ref * ones(N, 1); history_x zeros(n, sim_steps); history_u zeros(1, sim_steps); history_y zeros(1, sim_steps); for k 1:sim_steps xi [x; u_prev]; % 目标函数 H 2 * (G * Q * G R); f -2 * G * Q * (Y_ref - F * xi); % 约束 A_cons [tril(ones(M,M)); -tril(ones(M,M))]; b_cons [u_max - u_prev; -(u_min - u_prev)]; % 求解 options optimoptions(quadprog, Display, off, Algorithm, active-set); [dU_opt, ~, ~] quadprog(H, f, A_cons, b_cons, [], [], ... du_min*ones(M,1), du_max*ones(M,1), [], options); dU dU_opt(1); u u_prev dU; % 系统动态 x A * x B * u; y C * x; u_prev u; % 记录 history_x(:, k) x; history_u(k) u; history_y(k) y; end % 绘图 t (0:sim_steps-1) * Ts; figure; subplot(2,1,1); plot(t, history_y, LineWidth, 1.5); hold on; plot(t, ones(1,sim_steps), --); xlabel(时间/s); ylabel(输出); subplot(2,1,2); stairs(t, history_u, LineWidth, 1.5); xlabel(时间/s); ylabel(控制量);这段代码跑下来你应该能看到输出大约经过几次调整后稳定到 1控制量被限制在 [-1,1] 之间而且变化速率被限在 ±0.2。如果没有速率约束控制量可能会一步跳满这就是增量约束的价值。4. 关键体会三种公式切换时的坑与Matlab调试技巧4.1 预测矩阵维度不匹配是最常见的问题我见过太多人把增广状态搞对之后预测矩阵 F 和 G 的维度跟 C_tilde 相乘出错。你只要记得一句话F 的每一行对应一个预测步行的长度等于增广状态的维度G 的每一行也对应一个预测步列数等于控制增量个数通常就是 M。如果 F 是 N×(nm)G 是 N×M那 Y_pred 就是 N×1Y_ref 也必须是 N×1。检查维度时直接用 size 打印出来逐一对比盯着公式猜要快得多。4.2 控制时域 M 小于预测时域 N 的矩阵生成秘技如果 M N请注意 G 矩阵里只有前 M 列有值从第 M1 列开始应该全为零。上面的代码里内层 j 循环只到 M且 if j i 保证了当 i 较小时j 太大对应的项不存在。这里有个容易忽略的点当 i 很大、j 很小的时候A_tilde^(i-j) 的幂次很大数值上可能会出现近乎奇异的项这时应该考虑用递推而不是矩阵幂生成 G。对于一般的稳定系统N 不超过 20 问题不大如果 N 很大建议把生成 G 的循环改写成递推累加避免大幂次导致精度损失。4.3 quadprog 的算法选择与目标函数缩放Matlab 的 quadprog 默认算法在不同版本下可能不同有些老版本默认 active-set新版本默认 interior-point-convex。对于 MPC 这个规模不大、约束线性且上下界简单的 QP两者都能解。但注意 interior-point 算法可能会在约束比较紧时给出边界附近的解而 active-set 对初始点敏感。我通常的做法是显式指定 active-set因为 MPC 的 QP 一般只有几百个决策变量active-set 迭代次数少行为也容易预测。另外如果 H 的条件数太大因为 Q 和 R 相差几个数量级建议先对 Q 和 R 做归一化比如都除以输出量纲的平方否则求解器容易报“H 必须是正定”的警告。4.4 参考轨迹的预知与实时性问题上面代码里 Y_ref 全程是常数这只是最简化。如果参考轨迹是时变的需要在每个采样周期更新 Y_ref 向量。对于可以预知的参考轨迹比如规划的轨迹可以直接把未来 N 步的参考值放进去这就是“前馈式”MPC能显著改善跟踪滞后。如果你的应用只能拿到当前参考值那就让 Y_ref y_ref(k) * ones(N,1)也能用只是跟踪会有纯滞后。这在运动控制里要特别注意因为轨迹跟踪的滞后往往就是追踪误差的主要来源。4.5 三种公式在Matlab里的实际调试体验非增广公式在实现时最大的坑在于构造 Ψ 矩阵时容易把 u(k-1) 的传播项搞错。我曾经在调试时漏了 Γ_u 那一项结果系统响应奇慢后来发现是少加了一项 B*u(k-1) 的累加项。增广公式就完全不会出这种问题因为状态里带着 u(k-1)所有项都在预测矩阵里自动处理了。带扰动的公式调试难度主要在观测器如果 d 的估计滞后MPC 的补偿作用会来得慢这时需要检查扰动模型的收敛速度。4.6 给新手的Matlab调试建议我自己的习惯是三步走第一步先把未加约束的MPC闭环仿真跑通。只保留目标函数不加入任何不等式约束。如果这一步都不稳定那就是模型或预测矩阵错了此时去检查 A_tilde 和 B_tilde。第二步加入控制幅值约束。观察控制量是否触碰边界触碰边界时闭环是否仍稳定。如果触碰边界后发散多半是约束矩阵 A_cons 写错了尤其是符号问题。第三步加入速率约束和控制时域缩减。速率约束是增量形式的专属优势加上后观察响应变慢的程度再调 R。每加一层约束都单独验证千万不要一上来全加上否则出了问题根本分不清是预测错误还是约束冲突。5. 参数整定的经验笔记与常见问题速查5.1 预测时域N、控制时域M和采样周期Ts怎么配先说采样周期它必须是系统主要动力学的 1/5 ~ 1/10。比如系统开环时间常数是 1s那采样周期取 0.1 ~ 0.2s 比较稳妥。Ts 太大会丢失两采样点之间的动态太小会让每一步控制增量都很小导致跟踪迟钝。预测时域 N 的物理含义是“看多远的未来”至少要覆盖系统动态的主要时间跨度。一个简单经验N * Ts ≈ 3 ~ 5 倍系统时间常数。比如时间常数 1sTs0.1s那 N 至少取 30 到 50。取小了MPC 只看近处约束处理和稳定性都会变差。控制时域 M 决定了“自由度的数量”。M 越大控制越激进计算量也越大M 很小比如1控制量几乎只在第一拍决定后面全不更新响应会很保守。一般取 N 的 1/3 到 1/2 作为 M。5.2 Q和R的整定顺序跟 PID 调参一样先调 Q 让跟踪误差尽快收敛再用 R 抑制控制量的高频抖振。但要注意增量形式的 R 直接惩罚的是 Δu 的幅值R 太小控制增量会很大但因为有速率约束du_min/du_max兜底一般不会失控只是会让 u 频繁触及约束边界。R 太大系统响应变慢而且容易产生一个现象叫“提前减速”——MPC 预测到未来误差不大所以提前开始减小控制增量导致跟踪末端有明显的拖尾。遇到这种情况把 R 减半试试。5.3 约束冲突怎么办当输出约束、控制幅值约束和控制速率约束同时存在时可能出现“无可行解”。最简单的手段是扩大控制时域 M或者放宽速率约束。另一种更工程化的做法是将输出约束转成软约束引入松弛变量在目标函数里加一个惩罚项。Matlab里实现软约束比较简单你只需要把原问题里 A_cons 对应的输出约束改成 A_cons * dU b_cons slack同时在目标函数里加上 lambda * slack。这里 slack 是额外的决策变量H 矩阵对应位置为零f 向量对应位置加 lambda。这样无论约束是否冲突求解器都能给出一个解只是当约束被突破时slack 变大目标函数上升表现为“尽量满足约束但允许在极端情况下微量突破”。5.4 常见问题速查表现象可能原因排查方法闭环发散A_tilde 或 B_tilde 构造错误打印闭环极点检查增广矩阵的特征值仿真速度极慢每次循环里重复构造 F、G时不变系统 F、G 移到循环外输出稳态误差不为零模型无积分作用或没有扰动补偿改用带扰动增广公式或检查权重 Q控制量剧烈振荡R 太小或速率约束太松增大 R或收紧 du_min/du_maxquadprog 报错“H 不正定”R 为零或 Q 与 R 数量级悬殊给 R 加一个很小的正则项比如 1e-6参考轨迹跟踪滞后明显Y_ref 没有使用未来值用预知的参考轨迹填充未来 N 步控制量始终卡在下界约束写反了符号检查 A_cons 和 b_cons 的正负号5.5 一个调参案例我实际跑过的一个例子模拟项目X是一个温度控制系统采样周期0.5s系统时间常数约10s。最初取 N20M5Q1R0.1响应非常慢因为 N*Ts10s只覆盖了1倍时间常数MPC 看不到更远的未来。后来把 N 加到 40覆盖20sM 加到 10R 降到 0.02输出响应明显加快最终稳态误差为零。这就是“预测时域必须覆盖系统主要动态”的直观教训。还有一次我把速率约束设成 du_max0.5而实际系统每0.5s采样一次相当于控制量每秒变化不超过1这本身是合理的。但 R 设成了 0.001结果每步控制增量都碰到约束边界控制量一路滑到底系统超调接近30%。后来把 R 提到 0.05超调降到 8%。这说明速率约束和权重 R 相辅相成约束决定硬边界R 决定在边界内部的软取舍。6. 带扰动公式的Matlab扩展要点前面提到了第三种公式带扰动增广这里再补充一点代码层面的细节。扰动 d 通常不是直接可测的你需要用观测器估计。最简单的是一维扰动取扩展状态x_aug [x; u_prev; d]观测器设计可以用标准的扩展卡尔曼滤波器也可以用一个简单的龙伯格观测器前提是系统可观。Matlab里可以这样构造观测矩阵% 假设扰动模型: d(k1) d(k) % 观测方程对输出 y 的扰动输入是 C_d通常等于 [C, 0, 0; 0, 0, 1] 的形式 % 具体矩阵根据你的对象模型来定在MPC预测矩阵里F 和 G 需要对应扩展后的状态维度所以 F 的列为 n1_dim_u1_dim_dG 的列还是 M。扰动状态在预测中起到的作用是预测输出 名义模型输出 扰动在输出上的投影这样预测误差不再只是“当前误差”而是被扰动模型解释并持续到未来。实际调试中我建议先把扰动公式用在无扰动环境里验证确认闭环性能与经典增广完全一致再人为加入阶跃扰动观察MPC是否有补偿能力。如果发现补偿速度慢可以增大观测器增益但注意观测器增益过大会放大噪声这是一个需要权衡的点。7. 结束前的几个实操心得最后说点平时文档里不爱写的经验。第一做MPC仿真时千万别只盯着输出曲线和控制量曲线一定要把控制增量曲线也画出来。很多问题单看控制量觉得挺正常但控制增量已经高频抖动得不成样子这往往是 R 权重和约束不匹配的信号。我在做某图像处理Demo的辅助控制的时候就是靠画 Δu 曲线才发现一个采样周期里的数值振荡问题。第二代码模块化很重要。我说的是把“模型离散化”、“增广矩阵构造”、“预测矩阵生成”、“QP求解”、“仿真循环”分别写成独立的函数。初期调试用脚本没问题但一旦要换模型、换参数函数化能帮你节省大量时间。我自己的习惯是把预测矩阵生成写成 generatePredictionMatrices(A_tilde, B_tilde, C_tilde, N, M)输入输出一目了然后续切换公式时只需要替换这个函数内部实现。第三如果目标平台是嵌入式C代码Matlab里仿真用的 quadprog 通常不能直接移植。但你可以先用Matlab验证控制律然后把 QP 求解换成嵌入式实时求解器比如基于内点法或活动集法的轻量库。这个过程中增广公式的优势就体现出来了预测矩阵和约束矩阵的形式简单移植时不容易出错。相比之下非增广公式在移植时还得单独处理 u(k-1) 项容易埋bug。谁用谁知道。MPC没有想象中那么神秘但也没有某些论文写得那么轻巧。把输入增量这个点吃透不同状态空间公式对你来说就只是“怎么组织模型”的选择题而不是拦路虎。希望这篇能让你少走点弯路动手把代码跑起来之后你会发现MPC真正的乐趣在于“约束下求最优”的那个瞬间——一切看起来都在边界上跳舞但舞步是稳的。