MPC模型预测控制里优化变量到底是选 u 还是选 Δu这个问题我在刚入门时被绕了很久。后来在某伺服系统项目里做位置跟踪被稳态误差和教育了一回才真正意识到“输入增量”这四个字的含金量。现在回头看很多讲MPC的资料直接把增广状态空间公式甩出来却没有把“为什么要把输入变量改成增量”“不同公式之间到底差在哪”讲透。这篇文章就用Matlab代码实现的方式把我用输入增量构造不同的状态空间MPC公式的过程、推导步骤和踩坑记录整理出来主要面向正在自己手写MPC、不想只依赖现成工具箱的开发者。先说清楚范围我们只讨论离散线性状态空间模型目标是通过输入增量Δu作为优化变量实现几种不同的MPC公式包括直接在原状态空间上做增量优化、把输入扩成增广状态、以及输出增量型的velocity form。这三种公式最后都能化成同一个标准QP问题但预测矩阵、状态初始化和约束处理各有差异这正好对应标题里说的“不同的状态空间MPC公式”。我会把每一步的Matlab关键代码贴出来并给出我在仿真中踩过的坑和调参经验。1. 为什么MPC要用输入增量代替绝对输入先回到最基础的问题MPC在每个采样周期都要根据当前测量值用模型预测未来 N 步求解一个带约束的优化问题然后把得到的第一拍控制量送给被控对象。这里的优化变量如果取成绝对输入序列 u_k, u_{k1}, ..., u_{kN-1}就是我们常见的标准MPC形式取成增量序列 Δu_k, Δu_{k1}, ..., Δu_{kN-1}就变成了增量式MPC。从控制本质上看增量式 MPC 相当于在控制回路中自然引入了一个积分环节。原因很简单实际执行的控制量为 u_k u_{k-1} Δu_k控制器每次决策的是“在前一拍基础上加多少”。即使模型存在稳态误差或者有常值扰动只要闭环稳定积分作用会把输出慢慢拉回参考值。如果直接用绝对输入 u 作为优化变量相当于一个纯比例式的前馈控制方案遇到模型失配时很容易留下残余偏差。我用一个简单的类比普通MPC像你开车时只盯着当前位置和目标的距离打方向增量MPC则额外记住了方向盘上一拍在哪让你每一拍知道动了多少不会因为路面变化而越跑越偏。另外增量形式还能很自然地处理执行机构的速率约束。伺服电机、液压阀、舵机这些执行器通常都有自己的最大变化率限制也就是 |Δu| ≤ Δu_max。优化变量本身就是Δu那这个约束就直接写成变量的上下界处理起来非常方便。如果优化变量是绝对输入u你还需要单独再算相邻控制量之差麻烦不说还容易把预测方程搞得一团糟。但增量形式也有代价状态空间模型必须跟着变。原本的模型是 x_{k1} A x_k B u_k输入是 u现在输入变成了 Δu状态方程不会自动兼容。这时候就得对模型做改造可以是把上一拍输入 u_{k-1} 塞进状态向量也可以是用增量状态 Δx 和输出 y 组成一个新的增广系统。这就引出了多种所谓的“状态空间MPC公式”下面一章就逐个推导。2. 三种常见状态空间MPC公式推导为了方便说明先给定同一个离散线性定常系统x(k1)A x(k)B u(k)y(k)C x(k)其中 x 是 n_x 维状态u 是输入这里先讨论单输入情形y 是输出。A、B、C 可以是连续系统离散化得到也可以直接给离散模型。三种公式都假设当前时刻为 k已知测量状态 x_k、上一拍输入 u_{k-1}、上一拍状态 x_{k-1}某些公式需要。2.1 类型A原状态空间模型 控制增量作为优化变量第一种公式最简单状态空间模型保持原样不扩充状态变量只是把优化变量换成 Δu。虽然模型表面上没有变化但预测公式要额外处理 u 的累积关系。令预测时域 N_p N_c N优化变量为ΔU [Δu_k; Δu_{k1}; ...; Δu_{kN-1}]实际施加到系统的输入序列U [u_k; u_{k1}; ...; u_{kN-1}]根据 u_{kj} u_{kj-1} Δu_{kj}可以得到U M u_{k-1} L ΔU其中 M 是长度为 N 的全1列向量L 是 N×N 的下三角全1矩阵。比如 N3 时L [1 0 0; 1 1 0; 1 1 1]这一步是整个类型A的核心很多初学者在这里对不上矩阵维度。接着用标准预测公式Y F x_k G U其中 F 和 G 由 A、B、C 和预测时域决定第 j 行对应第 j 步预测输出。把 U 表达式代进去Y F x_k G (M u_{k-1} L ΔU) F x_k G M u_{k-1} G L ΔU这样预测输出就显式地表示成了优化变量 ΔU 的线性函数。代价函数取J (Y - Y_ref)^T Qbar (Y - Y_ref) ΔU^T Rbar ΔU这里 Qbar、Rbar 是对预测时域内的输出误差和增量输入加权的块对角矩阵。注意我没有把绝对输入 u 放进代价项因为既然已经用增量控制了一般不需要再单独惩罚 u如果你想限制控制力度可以在 QP 中额外加入约束或者再加一个 U 的惩罚项但这会让 Hessian 结构更复杂通常没必要。类型A的优点是状态维度没有增加QP 的变量数最少数值上比较轻量缺点是预测矩阵要额外地乘一个 L 矩阵而且 u_{k-1} 作为偏置项出现在预测公式里容易在代码实现时搞错下标。2.2 类型B输入增广状态空间模型第二种做法是干脆把上一拍输入 u_{k-1} 也写成状态变量组成一个新状态v_k [x_k; u_{k-1}]输入是 Δu_k。推导如下x_{k1} A x_k B u_k A x_k B (u_{k-1} Δu_k) A x_k B u_{k-1} B Δu_ku_k u_{k-1} Δu_k写成矩阵形式v_{k1} [[A, B], [0, 1]] v_k [[B], [1]] Δu_ky_k [C, 0] v_k这个增广模型非常整齐把输入累积关系直接嵌进了状态方程。预测输出可以写成Y F_b v_k G_b ΔU这里的 F_b、G_b 直接用标准预测矩阵构造函数生成不需要像类型A那样手工管理 L 和 M。代价函数形式和类型A一样。类型B的优势是模型统一后续扩展时变系统、加约束都比较方便而且不需要在预测公式里单独处理 u_{k-1} 的偏置因为 u_{k-1} 已经包含在 v_k 里。代价是状态维度从 n_x 变成 n_x1增广矩阵传递函数多了一个积分极点要求原系统可观测且增广系统可控否则会出现不可镇定模式后面仿真时要注意。2.3 类型C完全增量状态空间模型velocity form第三种公式更彻底不只把输入增量连状态增量 Δx 也纳入新的状态向量构成所谓 velocity form。定义w_k [Δx_k; y_k]其中 Δx_k x_k - x_{k-1}。推导过程稍微绕一点Δx_{k1} x_{k1} - x_k (A x_k B u_k) - (A x_{k-1} B u_{k-1}) A Δx_k B Δu_k这一步用到了 x_k x_{k-1} Δx_k 以及 u_k u_{k-1} Δu_k。再看输出y_{k1} C x_{k1} C (x_k Δx_{k1}) C x_k C(A Δx_k B Δu_k)因为 y_k C x_k所以y_{k1} y_k C A Δx_k C B Δu_k于是增广状态方程为w_{k1} [[A, 0], [C A, 1]] w_k [[B], [C B]] Δu_k输出直接从 w 中取第二项y_k [0, 1] w_k这里为了书写简单假设输出 y 是一维标量如果是多输出右下角的“1”要换成 I_n_y最后输出取对应块。这个模型里的输出本身成了状态积分作用最强对阶跃参考跟踪和常值扰动抑制天然无偏也是我在项目里最常用的形式。但类型C有两个实现上的坑。第一初始状态 w_k 需要 Δx_k x_k - x_{k-1}如果仿真第一拍把 x_{k-1} 错误设成 0起步就会产生一个很大的瞬时扰动脉冲第二如果 x 中某些分量不可直接测量需要设计状态观测器否则 Δx 的真实值无法获得。很多教材忽略了这个工程细节直到上手仿真才会被坑到。3. 三种公式如何统一写成标准QP不管是用上面哪种公式MPC最终都要在Matlab里用 quadprog 求解一个凸二次规划问题。通用形式是min 0.5 z^T H z f^T zs.t. A_eq z b_eqA_in z ≤ b_inlb ≤ z ≤ ub这里的 z 就是 ΔU 或与其等价的优化变量向量三种公式的差别只在于构造 H、f 和约束矩阵的方式不同。3.1 预测矩阵构造先写一个通用的预测矩阵构造函数避免每次重复贴一大串循环。函数根据 A、B、C 和预测时域 N返回两个矩阵 F 和 Gamma使得输出预测满足Y F x_k Gamma U_seq其中 U_seq 是绝对输入序列。代码如下function [F, Gamma] buildPrediction(A, B, C, N) ny size(C, 1); nx size(A, 1); F zeros(ny*N, nx); Gamma zeros(ny*N, N); for j 1:N F((j-1)*ny1:j*ny, :) C * A^j; for i 1:j Gamma((j-1)*ny1:j*ny, i) C * A^(j-i) * B; end end end注意这里的 F 中用的是 A^j 而不是 A^(j-1)因为预测的第一行对应 k1 时刻y_{k1} C(A x_k B u_k)。Gamma 的第 i 列对应第 i 个输入分量 u_{ki-1} 对输出预测的贡献所以第一行第一列是 C*B。我刚写这段代码时在这里习惯性地写了 A^(j-i-1)结果预测曲线像抽风一样乱跳后来一步一步手推才找出来。3.2 类型A的QP矩阵构造类型A中最终预测输出为Y F x_k Gamma * (M * u_last L * ΔU)记 G_A Gamma * Lbias_A F x_k Gamma * M * u_last则H_A 2 * (G_A * Qbar * G_A Rbar)f_A 2 * G_A * Qbar * (bias_A - Y_ref)约束包含四部分ΔU 上下界直接给 lb、ub绝对输入幅值约束u_min ≤ u_last L*ΔU ≤ u_max输出约束y_min ≤ bias_A G_A*ΔU ≤ y_max等式约束一般没有除非做终端约束或参考跟踪扩展。Matlab核心代码L tril(ones(N, N)); M ones(N, 1); G Gamma * L; bias F * xk Gamma * M * u_last; H 2 * (G * Qbar * G Rbar); f 2 * G * Qbar * (bias - ref); Ain [L; -L]; bin [u_max - u_last; u_last - u_min]; lb du_min * ones(N, 1); ub du_max * ones(N, 1); options optimoptions(quadprog, Display, off); [du_opt, ~, exitflag] quadprog(H, f, Ain, bin, [], [], lb, ub, [], options);这里我把输出约束省了避免矩阵维度把文章撑爆实际应用时输出约束也是类似写法GΔU ≤ y_max - bias以及 -GΔU ≤ bias - y_min。3.3 类型B和类型C的QP矩阵构造类型B的预测可以统一写成Y F_b * vk G_b * ΔU其中 F_b、G_b 用同一个 buildPrediction 函数生成传入增广矩阵 Ab、Bb、Cb 即可。类型C同理传入 Ac、Bc、Cc。于是H和f就是标准公式H 2 * (G_b * Qbar * G_b Rbar)f 2 * G_b * Qbar * (F_b * w0 - ref)约束中的输入幅值约束仍然是绝对输入表达式如果用类型B或Cu_{kj} u_last Σ_{i0}^{j} Δu_{ki}依然可以用 L 和 M 构造所以这个约束和类型A完全一样。增量约束也还是变量上下界。区别只在初始状态类型B是 [x; u_last]类型C是 [Δx; y]。我比较推荐先用类型C跑通整个仿真流程因为积分作用带来的无偏性会让你对MPC有很直观的“成了”的感觉。等类型C稳定了再回头对比类型B和A你会发现三类代码结构几乎一样只是模型矩阵和初始状态不同这就是标题里“使用输入增量实现了不同的状态空间MPC公式”的实际含义。4. Matlab仿真对比实现纸上谈兵差不多了接下来看具体代码。我以一个模拟的二阶线性对象为例连续传递函数取G(s) 1 / (s^2 0.8s 1)采样周期 Ts0.1s零阶保持器离散化。预测时域 N20控制时域 N也用20先不搞NcNp的复杂情况。权重 Q10R1输入幅值约束 |u|≤1增量约束 |Δu|≤0.5。4.1 系统离散化与参数设置s tf(s); Gc 1 / (s^2 0.8*s 1); Ts 0.1; Gd c2d(Gc, Ts, zoh); [A, B, C, D] ssdata(Gd); nx size(A, 1); N 20; Q 10; R 1; Qbar kron(eye(N), Q); Rbar kron(eye(N), R); u_min -1; u_max 1; du_min -0.5; du_max 0.5; ref 1.0; % 阶跃参考设定 T_sim 50; % 仿真步数注意这里的 Q 和 R 是单输出单输入的标量权重Qbar 和 Rbar 是对预测时域内每一步重复的块对角矩阵。如果系统是多输入多输出权重就要按维度扩展。4.2 仿真循环主体由于三种公式只差在模型和初始状态我先以类型C为例给主循环x zeros(nx, 1); x_prev zeros(nx, 1); u_prev 0; y C * x; u_log zeros(T_sim, 1); y_log zeros(T_sim, 1); du_log zeros(T_sim, 1); for k 1:T_sim % 类型C增广状态 Ac [A, zeros(nx, 1); C*A, 1]; Bc [B; C*B]; Cc [zeros(1, nx), 1]; w0 [x - x_prev; y]; [Fc, Gc] buildPrediction(Ac, Bc, Cc, N); bias_c Fc * w0; H 2 * (Gc * Qbar * Gc Rbar); f 2 * Gc * Qbar * (bias_c - ref); L tril(ones(N, N)); Ain [L; -L]; bin [u_max - u_prev; u_prev - u_min]; lb du_min * ones(N, 1); ub du_max * ones(N, 1); options optimoptions(quadprog, Display, off); [du_opt, ~, exitflag] quadprog(H, f, Ain, bin, [], [], lb, ub, [], options); if exitflag 0 warning(k%d: QP求解失败使用默认增量0, k); du_opt zeros(N, 1); end du du_opt(1); u u_prev du; % 更新被控对象状态可在这里注入扰动 d 0; x_new A * x B * u d; y_new C * x_new; u_log(k) u; y_log(k) y; du_log(k) du; x_prev x; x x_new; u_prev u; y y_new; end把 Ac、Bc、Cc 定义放在循环外面更好我这里放在循环内纯粹为了演示每步都能拿到最新的 w0。类型B只需要把增广矩阵改成Ab [A, B; zeros(1, nx), 1]; Bb [B; 1]; Cb [C, 0]; w0 [x; u_prev];类型A不需要增广矩阵直接按3.2节的公式构造 G 和 bias。三种代码跑完把 y_log 和 u_log 画出来对比你会发现一个很有意思的现象在模型完全准确、没有扰动的情况下三条响应曲线几乎重合一旦我在状态更新里注入一个常值扰动 d0.1类型A的输出会出现一个固定偏差而类型B和类型C凭借内部积分作用能把输出重新拉回参考值。这个差别就是“输入增量”的直观体现。4.3 三种公式的对比结果我实际仿真时的数据大致如下但不同对象和权重下会有差别关键看趋势对比项类型A原状态增量优化类型B输入增广状态类型C完全增量状态状态维度nxnx1nx1稳态误差无扰动000稳态误差有常值扰动有偏差基本为0基本为0对执行器速率约束处理需要L矩阵累加变量上下界直接变量上下界直接起步对初始状态敏感度中中高需x_prevQP变量数NNN预测矩阵构造复杂度中低低类型C在实际使用中通常是最稳的因为输出被做成状态后即使预测模型存在小幅误差积分作用也会持续修正不会让输出长期漂移。类型B在数学上等效于类型A加了一个积分状态但预测矩阵更整齐。类型A虽然代码量也不复杂但容易在累积矩阵上犯错而且遇到模型失配时天然吃亏。5. 常见问题与调参经验自己在Matlab里手写增量MPC时我踩过不少坑这里挑几个最典型的说。5.1 quadprog报错Hessian不是正定矩阵增量MPC的H矩阵原始形式是 GQbarG Rbar。如果G的列线性相关或者Rbar很小会导致H奇异。最简单的做法是在H上加一个很小的单位阵H H 1e-8 * eye(N);这样能保证quadprog有解。但不要加太大的值否则等于偷偷把增量权重调大响应会变慢。另外如果 Q0输出误差完全不惩罚H可能只由Rbar决定虽然正定但会失去跟踪能力这类低级错误也要排查。5.2 类型C起步跳变类型C的初始状态是 w0 [x - x_prev; y]。仿真第一拍如果我把 x_prev 设成了0而 x 是初始状态那 Δx 就是一个很大的非零向量。这个非零状态经过预测矩阵放大后会让第一拍增量 Δu 输出一个离谱的值。解决办法是如果系统初始确实静止在某个工作点直接令 x_prev x令 Δx0。很多人会犯这个错因为直觉上“初始状态为0”好像没问题但类型C的模型关心的是状态变化量而不是绝对量。5.3 增量约束与输入幅值约束冲突当 u_prev 已经在幅值上限附近而参考输入要求继续增大那么允许的 Δu 上限是 u_max - u_prev可能比 du_max 还小。我在代码中把绝对输入约束写成 L*ΔU ≤ u_max - u_last同时又把 lb/ub 设成 du_min/du_max这会同时生效。如果两个约束有冲突quadprog会报不可行。解决办法有两种一是把幅值约束和速率约束同时保留但把边界设置放宽二是用软约束给不可行留出余地。实际工程中执行器物理限制不能妥协所以通常会引入松弛变量。5.4 预测矩阵索引搞错前面说过 buildPrediction 的 Gamma 第一行第一列应该是 CB。很多人从书上的标准公式直接抄结果写成 CA^(j-i-1)B导致预测输出少了一拍。一个快速检查方法设 N1此时 F 应该是 CAGamma 应该是 C*B。如果输出不对大概率是索引问题。另外当 NcNp 时控制时域之后的 Δu 要补零L矩阵的列数小于行数很多初学者在这里维度报错。我一般建议先让 NcNp 跑通再去优化控制时域。5.5 权重Q和R的调参规律增量MPC中 R 直接约束的是 Δu 而不是 u所以 R 的物理意义和绝对输入MPC不一样。R 越大控制器越不愿意让输入变化太猛响应变慢但执行器寿命更长。Q 越大输出误差惩罚越重响应更快但容易出现振荡。我常用的试凑顺序是先把 R 设成一个和输出量级相近的值Q 设成 R 的5~10倍看阶跃响应如果超调太大增大 R如果稳态太慢减小 R。因为增量形式自带积分不用担心稳态误差所以调参压力比绝对输入MPC小很多。5.6 模型失配和扰动下的鲁棒性如果你用类型A发现稳态有偏差不要怀疑是代码问题这就是理论上的短板。类型B和类型C因为把积分环节放进了回路对模型失配的容忍度明显更高。但要注意积分作用也不是万能的如果模型误差导致闭环极点越过单位圆系统照样发散。因此实际项目中我会先用类型A做开环预测验证模型是否基本准确再用类型C做闭环控制这样排查问题更快。6. 个人体会与后续扩展方向写这篇文章的时候我又把三种公式的代码重新跑了一遍。最大的体会是不要被“不同的状态空间公式”这几个字吓到它们本质上是同一件事的三种姿态。类型A是在原模型上强行用增量变量类型B是把上一拍输入存进状态类型C则把状态增量也纳入了模型。选哪一种主要看你的系统状态是否容易获得、是否需要鲁棒抗扰动、以及你对预测矩阵整洁度的偏好。我个人在实际项目里首选类型C因为它最接近“把积分写进MPC”的思想后续加输出约束、软约束、参考轨迹滤波都顺手。如果还想继续深入建议往三个方向扩展一是把单输入单输出改成多输入多输出这时 Qbar 和 Rbar 要用 kron 正确处理维度输出约束也要按通道分别给出二是做非线性系统的增量式线性时变MPC在每个采样点重新线性化得到新的A、B矩阵再用类型B或C的框架滚动求解三是加终端代价或者终端约束理论上可以保证闭环稳定性但这需要另外设计终端权重矩阵工程实现会复杂不少。最后再分享一个调试小技巧写增量MPC时先不执行 u u_prev du而是先把 du 和 u_last 分别画出来观察是否出现剧烈振荡。如果 du 在拼命挣扎但输出没有改善大概率是预测方向反了或者状态初始化有误。这个检查习惯帮我省下了很多排查时间。