最近和几位做过程控制的朋友聊到一个很有意思的问题模型预测控制Model Predictive Control, MPC里为什么很多人不用u做决策变量非要绕一圈用输入增量Δu更让人头大的是不同教材和论文里给出的状态空间MPC公式长得还不一样有的用状态扩展有的用输出误差增广到底哪个能直接搬到Matlab里跑我花了两周时间把一套基于输入增量的状态空间MPC从公式推导到Matlab代码完整跑通顺便把两种主流公式做了横向对比。这篇把中间的所有细节捋一遍包括推导过程、矩阵构造、quadprog求解还有我踩过的几个坑给正在入坑MPC的同学做个参考。1. 问题引入为什么MPC要用输入增量公式1.1 直接控制输入会遇到什么尴尬如果你刚接触MPC最自然的想法是把未来N步的控制输入序列u(k)、u(k1)...u(kN-1)直接当作决策变量目标函数里惩罚输出误差然后求解。这种做法在理论上是成立的很多简化教材也这么讲但一到了实际工程就露馅了。举个例子现场有一个调节阀输入u是阀门开度范围0到100%。如果只把u当作决策变量MPC的优化器可能给出这样一组指令第1步开度从40%跳到90%第2步又从90%跳到10%。这在数字仿真里看着没问题但阀门的执行机构受不了连杆和阀芯会很快磨损甚至会引发管道压力波动。你想限制阀门每分钟最多只能变化5%这时候如果决策变量是u你只能去构造一个很别扭的约束|u(k)-u(k-1)| ≤ 5。虽然数学上也能表达但每次新增一个变化率约束都要重新拆矩阵非常麻烦。而如果把Δu u(k) - u(k-1)作为决策变量那么变化率不能超过5直接变成Δu的上下界约束-5 ≤ Δu ≤ 5。幅值约束和变化率约束在同一个框架下被统一处理代价函数里也可以直接对Δu加权相当于对控制动作的猛烈程度做惩罚这是增量式MPC第一个直观优势。1.2 增量式天然带了一个积分器第二个优点更隐蔽但更重要增量式MPC在控制回路里悄悄插入了一个离散积分器。为什么因为最终作用于被控对象的u(k)等于上一次的u(k-1)加上当前决策的增量Δu(k)也就是u(k) u(k-1) Δu(k)从频域看这相当于在控制器输出端串联了一个离散积分环节。这意味着即使被控对象本身是0型系统无积分环节整个开环回路也包含了一个积分因子闭环系统对阶跃型模型失配和外部常值扰动具有天然的抑制能力稳态误差可以减到零。这在直接以u为决策变量的MPC里往往需要额外引入扰动估计或输出误差积分项才能做到。我的实际感受是增量式MPC解决的不只是约束怎么写的问题它顺带把无静差这个控制工程师非常关心的问题一并解决了。这也是很多工业MPC软件默认使用增量式的原因。1.3 增量式与普通状态空间公式的定位差异既然有了增量决策变量原来的状态空间模型就不能直接用了。你需要把增量关系嵌入到模型中构造出新的状态空间描述。不同教材给出了不同的构造方法本质上都是在状态里多放一些历史信息。这里先记住一个概念所谓输入增量MPC核心是模型中的控制变量从u变成了Δu而状态变量经过了重新组合可能是[x; u(k-1)]也可能是[Δx; e]。有了这个基础后面所有矩阵推导都不会乱。2. 理论推导从原始状态空间到增量状态空间2.1 原始模型的基准坐标先写清楚最普通的线性时不变离散状态空间模型x(k1) A x(k) B u(k) y(k) C x(k) D u(k)其中x维度n×1u维度m×1y维度p×1。注意D是直接馈通矩阵物理系统中大多数情况下D0但推导时我先保留它后面会讨论D≠0带来的麻烦。为了构造增量模型我们把新的决策变量定义为Δu(k)u(k)-u(k-1)。直接拿这个Δu去替换原方程里的u是不行的因为原方程里没有历史输入项。所以第一个任务是把u的历史值变成状态的一部分这样表达式里才能出现Δu。2.2 状态扩展法把u(k-1)塞进状态取新的状态向量ξ(k) [x(k); u(k-1)]。也就是说把上一拍的控制量作为当前状态的一部分。那么下一拍的状态由两部分组成x(k1) A x(k) B [u(k-1)Δu(k)] A x(k) B u(k-1) B Δu(k) u(k) u(k-1) Δu(k)写成增广矩阵形式ξ(k1) [A B; 0 I] ξ(k) [B; I] Δu(k)输出方程变为y(k) C x(k) D [u(k-1)Δu(k)] [C D] ξ(k) D Δu(k)如果假设D0输出方程会干净很多变成y(k) [C 0] ξ(k)。这也是大多数增量MPC实现里的默认假设。这里的增广矩阵A_e [A B; 0 I]B_e [B; I]C_e [C 0]D0时。维度从原来的n变成nm增加的维度正好是输入通道数。这样做的好处是模型结构清晰可以处理非恒定参考轨迹、可以方便地加入输入幅值约束坏处是状态里多了一堆u的历史值噪声环境下估计这些状态需要额外小心。2.3 误差增广法差分加跟踪误差另一种思路是从系统差分出发得到一组以Δx和输出误差e为状态的模型。先定义两个差分变量Δx(k) x(k) - x(k-1) Δu(k) u(k) - u(k-1)对原始状态方程做差分Δx(k1) x(k1)-x(k) A[x(k)-x(k-1)] B[u(k)-u(k-1)] A Δx(k) B Δu(k)输出误差定义为e(k) y(k) - r其中r是恒定参考目标。如果D0那么e(k1) y(k1)-r y(k)-r [y(k1)-y(k)] e(k) C[x(k1)-x(k)] e(k) C[AΔx(k)BΔu(k)]写成增广矩阵[Δx(k1); e(k1)] [A 0; CA I] [Δx(k); e(k)] [B; CB] Δu(k)输出量就是e(k)即e(k) [0 I] [Δx(k); e(k)]注意这个模型中状态维度是np状态差分加输出误差不再需要存u(k-1)。它的含义非常直观状态的一部分是状态的增量另一部分是输出对参考的误差目标函数里直接惩罚e控制器天然就是一个积分型控制器。这种形式在推导预测方程的时候也比较舒服因为优化问题要跟踪的参考r不直接出现在输出方程里误差e本身已经是跟踪目标。不过它也有一个限制要求参考轨迹恒定或者至少能预先知道参考差分Δr否则e(k1)的表达式里要额外减去Δr(k1)模型就变复杂了。2.4 两种公式等价吗结论是在相同工况下两种公式最终闭环行为可以实现一致因为它们描述的是同一个物理系统只是状态坐标不同。状态扩展法把历史输入当状态通用性强误差增广法把输出误差当状态直观性更强更适合恒定参考跟踪。实际项目里我一般优先用状态扩展法因为遇到变参考、扰动和D不为零的情况不用改模型结构。只有在写学术报告时才特意用误差增广法展示无静差特性。3. 不同MPC公式的具体形式与实现要点3.1 基于状态扩展式的预测方程无论采用哪种增广模型MPC的最后一步都要把未来Np步的预测输出表示成当前状态和未来Δu序列的线性函数。以状态扩展式为例假设D0模型为ξ(k1) A_e ξ(k) B_e Δu(k) y(k) C_e ξ(k)令Np为预测时域Nc为控制时域且Nc≤Np未来控制增量为ΔU [Δu(k); Δu(k1); ...; Δu(kNc-1)] 维度Nc*m×1迭代展开y(k1) C_e A_e ξ(k) C_e B_e Δu(k) y(k2) C_e A_e^2 ξ(k) C_e A_e B_e Δu(k) C_e B_e Δu(k1) ... y(kNp) C_e A_e^{Np} ξ(k) C_e A_e^{Np-1}B_e Δu(k) ... C_e A_e^{Np-Nc}B_e Δu(kNc-1)把输出堆成列向量Y [y(k1); y(k2); ...; y(kNp)]则Y Γ ξ(k) Θ ΔU其中Γ是Npp×n的矩阵由C_e A_e^i堆叠而成Θ是Npp×Nc*m的下三角块Toeplitz矩阵。这一步是代码实现的核心。3.2 基于误差增广式的预测方程误差增广式的推导过程完全一样只不过状态变量替换成xa [Δx; e]输出是e参考项不出现在预测方程里。最后的预测式形式相同E Γ_a xa(k) Θ_a ΔU其中E是预测误差向量目标函数直接写成J E^T Qbar E ΔU^T Rbar ΔU。由于e已经是误差不需要像状态扩展式那样在目标函数里减去参考轨迹r代码实现少一个变量。3.3 输入幅值约束的矩阵化这是实战里最需要小心的地方。增量MPC的决策变量是ΔU但现场约束往往给的是u的上下限需要把两者关系转成矩阵不等式。因为u(ki) u(k-1) Σ_{j0}^{i} Δu(kj)把所有Nc步的控制量写出来U I_u u(k-1) L ΔU其中I_u是Nc*m×1的全1? 准确说每个控制通道有一个初始值L是一个由0和1组成的下三角块矩阵。例如单输入Nc3时U [u(k); u(k1); u(k2)] [u(k-1); u(k-1); u(k-1)] [1 0 0; 1 1 0; 1 1 1] ΔU所以幅值约束umin ≤ U ≤ umax变成两组线性不等式L ΔU ≤ umax - u(k-1)·ones -L ΔU ≤ -umin u(k-1)·ones直接把这些不等式堆进quadprog的不等式约束里。增量约束本身就是上下界lb/ub不需要额外处理。状态约束比如输出上限则需要通过预测方程YΓξΘΔU再进行线性变换实际操作时把约束写成ΔU的线性不等式即可。3.4 终端项和稳定性严格来说有限时域MPC不加终端约束或终端惩罚只能靠足够大的Np来近似保证稳定性不能给出确定性的稳定性证明。工程上常见做法是在代价函数叠加一个终端惩罚项E(kNp)^T P E(kNp)其中P可以从离散李雅普诺夫方程或黎卡提方程取。对增量MPC我建议至少加一个小的终端惩罚能明显抑制预测末端的振荡在数值上也更好看。4. Matlab代码实现详解4.1 构建预测矩阵的自定义函数我写的代码不使用MPC工具箱纯手算矩阵再调quadprog这样每一步都能看得见。先给一个通用的预测矩阵构建函数function [Gamma, Theta] buildPrediction(Ae, Be, Ce, Np, Nc) n size(Ae, 1); p size(Ce, 1); m size(Be, 2); % Gamma: Np*p x n Gamma zeros(Np*p, n); Ai Ae; for i 1:Np Gamma((i-1)*p1:i*p, :) Ce * Ai; Ai Ai * Ae; end % Theta: Np*p x Nc*m Theta zeros(Np*p, Nc*m); for j 1:Nc % 第j个控制增量对第i步预测的影响 Mi Be; % A_e^{i-j} * B_e, 初始为 B_e for i j:Np Theta((i-1)*p1:i*p, (j-1)*m1:j*m) Ce * Mi; Mi Ae * Mi; end end end这个函数假设可以一次性把Np步预测和Nc步控制都算出来没有考虑Nc后续预测中Δu保持为0的情况。如果你设了控制时域Nc小于预测时域Np上面Theta矩阵中第jNc的列保持为0表示Nc步之后增量不再变化。使用前先在命令窗口打印size确认和预期一致。4.2 二次规划问题装配有了Gamma和Theta目标函数可以变成标准的二次型。首先定义输出权重矩阵Q和增量权重矩阵R。如果p1Qbar kron(eye(Np), Q)Rbar kron(eye(Nc), R)。多输出时Q为p×p矩阵R为m×m矩阵构建方法类似。对于状态扩展式参考轨迹是ref向量长度为Npp预测误差E_ref Gammaz0 - ref其中z0是当前扩展状态。目标函数写成J E_ref^T Qbar E_ref ΔU^T Rbar ΔU展开得到J ΔU^T (Theta^T Qbar Theta Rbar) ΔU 2 E_ref^T Qbar Theta ΔU 常数项quadprog要求目标函数形式为0.5*x H x f x因此H_qp 2 * (Theta * Qbar * Theta Rbar); f_qp 2 * Theta * Qbar * E_ref;注意E_ref Gamma*z0 - ref所以f_qp是列向量。约束部分先用3.3的公式构建Aineq和bineq再设定lb和ub为每个通道的Δu下限/上限。如果没有增量约束用空数组[]。完整在线求解核心代码% 设置系统模型示例2x2系统 A [1.1 0.1; 0 0.8]; B [0.1 0; 0.2 0.1]; C [1 0; 0 1]; D zeros(2,2); % 离散模型已经是离散的直接采样时间Ts1 % 状态扩展矩阵 Ae [A, B; zeros(size(B)), eye(size(B,2))]; Be [B; eye(size(B,2))]; Ce [C, zeros(size(C,1), size(B,2))]; Np 10; Nc 4; [Gamma, Theta] buildPrediction(Ae, Be, Ce, Np, Nc); Q eye(2); R 0.1*eye(2); Qbar kron(eye(Np), Q); Rbar kron(eye(Nc), R); H_qp 2 * (Theta * Qbar * Theta Rbar); f_qp_base 2 * Theta * Qbar; % 状态初值和历史输入 x0 [0; 0]; u_prev [0; 0]; z0 [x0; u_prev]; ref [1; 0.5]; % 恒定参考 ref_seq repmat(ref, Np, 1); options optimoptions(quadprog, Display, off); for k 1:50 E_ref Gamma*z0 - ref_seq; f_qp f_qp_base * E_ref; % 输入幅值约束每个通道下限0上限[1.2;0.8] L tril(ones(Nc)); L_block kron(L, eye(2)); Aineq [L_block; -L_block]; bineq [1.2*ones(2*Nc,1) - kron(ones(Nc,1), u_prev); 0*ones(2*Nc,1) kron(ones(Nc,1), u_prev)]; lb -0.1*ones(2*Nc,1); ub 0.1*ones(2*Nc,1); dU quadprog(H_qp, f_qp, Aineq, bineq, [], [], lb, ub, [], options); if isempty(dU) error(优化器无解); end du dU(1:2); u u_prev du; x A*x0 B*u; % 实际系统 y C*x; % 更新状态 x0 x; u_prev u; z0 [x0; u_prev]; % 记录数据 end这里有个小细节约束里的下限用了0*ones是为了把u的下限0和上限[1.2;0.8]统一写成bineq。实际系统如果有差异按3.3的公式替换就好。4.3 闭环仿真结果观察用上面的2输入2输出示例系统跑完50步参考设为[1;0.5]输入约束上限为[1.2;0.8]增量约束在±0.1内可以看到两个输出均在15步内收敛到参考稳态误差为零因为输入增量模型本身带积分结构。把增量约束从±0.1放宽到±0.3收敛速度明显加快但控制动作更激进。这是增量权重的直观体现。可以做一个简单的对比表格供调试时选择参数| Np | Nc | 输出峰值超调 | 调节时间步 | 最大|Δu| | 计算耗时 | |----|----|--------------|----------------|------------|----------| | 10 | 3 | 8% | 14 | 0.35 | 12ms | | 20 | 5 | 4% | 9 | 0.28 | 20ms | | 30 | 5 | 3% | 7 | 0.26 | 31ms |以上数据来自我测试用的一个特定模型具体数值会因为对象不同而变但趋势是普遍的Np越长预测越远超调越小但计算量线性增加Nc太长会导致优化变量变多且容易出现病态。4.4 误差增广式代码差异如果换成误差增广式预测矩阵的构建过程完全一样只需要把A_e、B_e、C_e替换成A_inc [A, zeros(2, p); CA, eye(p)]; B_inc [B; CB]; C_inc [zeros(p, n?)?];具体代码中要保证维度正确。例如对上面2输入2输出系统p2n2则A_inc维度4×4B_inc 4×2C_inc 2×4。初始状态xa0 [0;0; -ref_target; 0? 初始误差e0 y0 - ref -ref]。然后把参考序列ref_seq从目标函数中去掉直接用预测误差向量。个人经验是一旦状态扩展式跑通误差增广式只是换几个矩阵核心逻辑不变。5. 常见问题与排查技巧5.1 预测矩阵维度不匹配这是新手最容易卡住的地方。我见过最常见的错误是Theta的行数写成Np而不是Npp或者列数写成Nc而不是Ncm。排查方法很简单在构建预测矩阵后立刻打印size(Gamma)和size(Theta)多输出系统尤其要检查。如果Gamma第一维不是Np*p说明Ce的行数或者堆叠语句写错了。更隐蔽的问题是控制时域Nc小于预测时域Np时Theta的后Np-Nc列仍然是零矩阵这部分是对未来增量不再变化的体现。如果不想要这个假设需要把Nc设为等于Np或者单独处理。5.2 quadprog返回无解增量MPC出现不可行问题非常常见典型情况是输入幅值约束和增量约束互相矛盾。比如u当前是1.0上限是1.0而增量下界是-0.05但优化器在某个预测步发现要把u增大到1.05才能跟踪参考结果被幅值约束卡住约束集合变成空集。这时候优先检查bineq是否包含了u(k-1)的偏移项kron(ones(Nc,1), u_prev)如果漏掉约束矩阵会整体偏掉一个初始值。如果真的找不到原因可以给约束加软变量。对输出约束和输入约束都可以引入松弛变量ε目标函数增加权重ρ*ε^2这样即使约束冲突优化器也能给出一个尽量满足的解而不是直接返回空。工程上这种做法几乎必须要有因为现场工况千变万化硬约束太容易导致控制器罢工。5.3 直接馈通矩阵D带来的输出方程额外项当原系统D≠0时状态扩展式的输出方程多了一个DΔu项预测公式里Theta会多出一块。很多教材默认D0导致读者在实际系统遇到直接馈通时一仿真就发现预测输出和真实输出对不上。我的建议是如果原始模型D确实不为零先不要直接套模型尽量通过合理建模把D消掉。比如传感器和执行器之间如果有纯延迟或滤波器一般可以把D置零并调整B/C矩阵。如果无法消除需要在预测方程的Theta里把D的影响叠加上——具体做法是定义一个新的C_tilde [C D]并且要特别小心第一行预测输出因为它会受到当前Δu的直接影响。这时Theta矩阵的第一行对应C_tilde*Be? 注意由于D直接乘Δu表达式与Be块会有交叉实际推导时按原始状态方程迭代不要直接沿用D0的Theta公式。5.4 参考轨迹非恒定状态扩展式可以直接处理随时间变化的参考轨迹只要把未来的ref序列准备好塞进E_ref Gamma*z0 - ref_seq即可。而误差增广式需要参考差分项Δr模型会变成e(k1) e(k) C A Δx(k) C B Δu(k) - Δr(k1)如果你不想改模型可以对参考信号做前馈补偿或者直接把参考序列设为零均值偏差。更省事的做法是非恒定参考路径一律用状态扩展式恒定参考用误差增广式各取所长。5.5 数值病态与权重太小时的处理增量权重R若设置太小H矩阵可能变成接近奇异quadprog求解时会出现警告甚至乱七八糟的结果。我的处理方式是先给R一个合理的初值比如0.01~0.1然后逐步减小观察控制动作的变化。如果发现输出有高频抖振多半是R太小而不是控制器频带不够。也可以在H上加一个正则化项R_eps 1e-8 * eye(Nc*m)相当于让每个控制增量付出极小代价。这个方法在矩阵阶数高、Q权重很大时非常有效。5.6 启动阶段的u(k-1)怎么定第一次进入MPC时u(k-1)是未知的。最稳妥的办法是把上一拍实际发给执行机构的指令值读回来。如果系统刚上电可以用一个安全的默认值或者直接令u(k-1)0同时把增量约束的初始界放宽让控制器在第一个周期内快速逼近合理值。6. 个人体会与扩展建议这套代码前前后后改了三天最深的体会是两种增量状态空间公式不是谁比谁高级而是适用场景不同。状态扩展式是万能钥匙什么工况都能套误差增广式是手术刀在恒定参考、无扰动的情况下公式简洁状态含义清楚调试起来特别舒服。我建议新入门的同学先用状态扩展式把整个MPC闭环跑通彻底理解Gamma和Theta怎么来再回头对比误差增广式你会发现很多概念一下子就能串起来。最后分享一个小技巧在Matlab里做MPC仿真时不要一上来就上四阶系统先用一个一阶惯性加纯延迟的简单模型把预测矩阵、约束映射、quadprog求解全部验证正确再扩展到多变量系统。这样排查问题的时间可以节省一大半。后续如果你想做非线性MPC或者分布式MPC这套增量状态空间的矩阵构造逻辑仍然是核心的底层框架不会白学的。