1. 为什么输入增量会成为MPC工程落地的分水岭刚接触MPC模型预测控制时很多人会直接照抄教材里的状态空间公式u(k) -K * x(k)或者在线优化时把控制量u本身直接当作决策变量丢进二次规划里。这种做法在仿真里跑双积分器、倒立摆完全没有问题但一旦拿去接真实被控对象十有八九会出问题执行器抖得像筛糠、稳态精度达不到、甚至系统直接发散。问题出在哪因为真实控制器输出的是绝对控制量而绝大多数执行机构和被控对象真正响应的是控制量的变化。电动阀门关心的是开度增量而不是开度本身伺服电机驱动器接受的是速度增量命令而不是速度绝对位置锅炉燃料阀在工况切换时怕的是燃料量突变而不是当前燃料量偏大。我最早挨过一刀的场景是温控系统采样周期30秒MPC算出来u(k)62.7%开度结果执行机构咔哒一声直接猛转温度超调了4度。后来把问题想透了——控制器本质上是在给每一个采样时刻下发一个新的绝对位置而机构响应的是这个位置的跳变。跳变越大冲击越大超调越严重。这就是标题里“输入增量”这个说法的工程来源。所谓输入增量Incremental Input本质上是把决策变量从u(k)换成Δu(k) u(k) - u(k-1)让优化器去计算“下一步该相对当前多开多少、少开多少”而不是直接计算“下一步开到哪里”。这个看似简单的变量替换在公式结构和求解难度上带来的是完全不同的两套MPC公式体系。Matlab里做这件事最直观的路线就是改写成增广状态空间形式把输入增量当成新系统的控制量然后再套标准的状态空间MPC求解框架。这篇文章就围绕这件事展开输入增量如何进入状态空间表达式、两种常见的增广公式写法、Matlab实现里的坑以及我怎么在仿真里验证这种改写带来的实际收益。2. 状态空间MPC里“控制量”和“控制增量”的本质差异2.1 原始状态空间MPC的表达形式离散状态空间模型的标准写法是x(k1) A·x(k) B·u(k) y(k) C·x(k) D·u(k)标准MPC的做法是在当前时刻k预测未来Np步的系统输出把未来Nc步的控制量u(ki)当作决策变量最小化一个形如J Σ ||y(ki) - r(ki)||²_Q Σ ||u(ki)||²_R的目标函数。约束可以加在u幅值上也可以加在输出y幅值上。这个框架在数学上是干净的但工程上有一个被掩盖的问题优化变量是绝对量u优化器不知道上一时刻的u(k-1)是什么。也就是说它计算出的最优解可能让u(k)比u(k-1)高出一大截——而这个一大截正是执行机构冲击的来源。你可以在约束里加|u(k)-u(k-1)| ≤ Δu_max来限制但这等于额外引入了一组约束不等式而且这种“软限制”并不能让优化器主动去利用增量信息做前瞻性的调节。2.2 输入增量作为决策变量的“状态增广”思路把增量引入状态空间MPC最通用、最经典的技巧是状态增广。具体来说把当前时刻的控制量u(k)也纳入状态向量构造新的状态变量ξ(k) [ x(k)^T, u(k-1)^T ]^T这样一来新系统的“控制输入”就变成了Δu(k) u(k) - u(k-1)新系统的状态方程可以推导为ξ(k1) [ x(k1) ] [ A B ] [ x(k) ] [ B ] Δu(k) [ u(k) ] [ 0 I ] [ u(k-1) ] [ I ]也就是说ξ(k1) A_tilde · ξ(k) B_tilde · Δu(k) y(k) C_tilde · ξ(k) D_tilde · Δu(k)其中A_tilde [ A B ] [ 0 I ] B_tilde [ B ] [ I ] C_tilde [ C D ] D_tilde D注意一个细节即使原始系统的D矩阵为0绝大多数物理系统直接传递项为0增广系统里y(k)仍然通过C_tilde中那个历史控制项[C D]·[x(k); u(k-1)]受到u(k-1)的影响也就是说输出预测天然包含了过去控制量的余效。这正是增量式MPC“记住历史位置”的数学体现。2.3 数学上的分水岭从控制幅值转为控制变化率两种公式的核心区别可以用一句话概括位置式MPC优化的是u(k)的绝对幅值目标函数里压制的是“控制量大小”。增量式MPC优化的是Δu(k)的变化率目标函数里压制的是“控制量变化速度”。在目标函数中增量式MPC的二次项变成了J Σ ||y(ki) - r(ki)||²_Q Σ ||Δu(ki)||²_R这里的R直接约束的就是执行机构的动作剧烈程度。工程上修改R的物理含义变得极其清晰——R越大Δu被压得越小阀门动作越柔和。位置式MPC里R压制的是绝对位置调参时很难直观理解“R0.1到底是让阀门开度小还是让阀门动作慢”而增量式MPC里R0.1直接告诉你“每一步阀门开度变化最多被惩罚到0.1的平方量级”。3. 两种增量式MPC公式的推导与对比从Matlab代码看差异3.1 方案一直接增广状态向量状态增广法我刚入行时在Matlab里是这样实现的。假设原始系统是二阶双积分器A [1 0.1; 0 1]; B [0; 0.1]; % 采样周期0.1s C [1 0]; D 0; % 状态增广 [nx, nu] size(B); A_tilde [A, B; zeros(nu, nx), eye(nu)]; B_tilde [B; eye(nu)]; C_tilde [C, D]; D_tilde D; % 新状态维度 nx_tilde size(A_tilde, 1);这个方案的好处是结构非常直白代码可读性高。预测矩阵F和Phi直接用增广后的A_tilde、B_tilde构建推导出来的每一步预测输出都可以写成一个标准的线性表达式。这个方案适合在控制对象维度比较低比如SISO系统、2-3阶系统时快速验证。3.2 方案二构建增量形式的扩展模型速度形式法另一种在学术文献和工业软件中更常见的写法是构造扩展状态变量将原状态x(k)、控制量u(k-1)、以及系统输出y(k)一起纳入状态然后在等式两边同时做差分Δx(k1) A·Δx(k) B·Δu(k) Δy(k1) C·Δx(k1) D·Δu(k1)再把y(k)作为增广状态的一部分得到[ Δx(k1) ] [ A 0 ][ Δx(k) ] [ B ] Δu(k) [ y(k1) ] [ CA 1 ][ y(k) ] [ CB ]这就是很多MPC教材里讲的“速度形式模型”Velocity-form Model它的关键特点是不再包含绝对状态x(k)而是全部用增量Δx(k)和输出y(k)来描述。这个写法在无静差跟踪问题上特别有效因为它天然嵌入了一个积分环节相当于在模型内部含有一个隐式积分器。Matlab代码实现上这种方案会稍显绕A_vel [A, zeros(nx, ny); C*A, eye(ny)]; B_vel [B; C*B]; C_vel [zeros(ny, nx), eye(ny)]; D_vel D;ny是输出维度。这个增广结构的核心是用y(k)替代了状态里的一部分预测输出直接用增广状态中的输出分量提取不需要额外乘C矩阵。3.3 方案对比什么时候用哪种对比维度方案一直接增广方案二速度形式状态物理含义原状态历史控制量状态增量历史输出稳态无静差需额外加积分补偿天然含积分作用代码可读性直观、好debug较绕、容易混淆维度输出预测提取需乘C_tilde直接取状态分量适合场景快速验证、教学工程跟踪、含扰动场合我在自己项目里通常先用方案一快速验证MPC核心逻辑是否正确确认无误后再切换到方案二去跑正式工况。原因很简单方案一的代码每一步变量名都很直观一旦运行结果不对用disp打印中间矩阵就能排查方案二一旦维度搞错错误信息会藏在增广矩阵内部排查成本高好几倍。4. Matlab实现一步步拆解从预测矩阵到quadprog求解4.1 预测矩阵的构建逻辑不管哪种增广方式MPC进入在线求解之前的核心工作都是构建预测模型。以方案一为例将增广后系统A_tilde、B_tilde代入标准MPC预测表达式未来Np步输出可以写成Y F·ξ(k) Φ·ΔU其中F [ C_tilde·A_tilde ] [ C_tilde·A_tilde² ] [ ... ] [ C_tilde·A_tilde^Np ] Φ [ C_tilde·B_tilde 0 ... 0 ] [ C_tilde·A_tilde·B_tilde C_tilde·B_tilde ... 0 ] [ ... ] [ C_tilde·A_tilde^(Np-1)·B_tilde ... C_tilde·B_tilde ]Matlab里我习惯用循环构建不推荐符号推导Np 20; % 预测时域 Nc 5; % 控制时域 F zeros(Np*ny, nx_tilde); Phi zeros(Np*ny, Nc*nu); % 构建F矩阵 for i 1:Np F((i-1)*ny1:i*ny, :) C_tilde * (A_tilde^i); end % 构建Phi矩阵 for i 1:Np for j 1:min(i, Nc) Phi((i-1)*ny1:i*ny, (j-1)*nu1:j*nu) ... C_tilde * (A_tilde^(i-j)) * B_tilde; end end注意Phi矩阵的列数是Nc*nu不是Np*nu——超过控制时域之后的输入增量保持为0这是“控制时域缩短”的基本思想。我第一次自己写的时候在这里踩过坑导致矩阵维度不匹配报错。4.2 二次规划目标函数与约束的拼装预测输出Y与参考轨迹R的偏差可以写成E R - F·ξ(k)目标函数J (R - Y)^T · Q̄ · (R - Y) ΔU^T · R̄ · ΔU展开后二次项矩阵H和一次项向量f为H Φ^T · Q̄ · Φ R̄ f -Φ^T · Q̄ · E其中Q̄是Np*ny × Np*ny的输出权重矩阵R̄是Nc*nu × Nc*nu的输入增量权重矩阵。Matlab里直接用kron构建Q_bar kron(eye(Np), diag([1 1])); % 输出权重 R_bar kron(eye(Nc), diag([0.1 0.1])); % 增量权重 H Phi * Q_bar * Phi R_bar; f -Phi * Q_bar * (Rs - F * xi_k);约束方面增量式MPC可以同时处理三类约束% Δu幅值约束 lb -0.5 * ones(Nc*nu, 1); ub 0.5 * ones(Nc*nu, 1); % u幅值约束通过累积和矩阵处理 A_cons tril(ones(Nc*nu)); b_lb -u_min u_prev_rep; % 累积下限 b_ub u_max - u_prev_rep; % 累积上限这里的u_prev_rep是把u(k-1)重复扩展到Nc步的向量。对于SISO系统nu1来说tril(ones(Nc))矩阵把Δu从1到Nc步累加就等于每一步的绝对控制量。这个约束矩阵的处理是增量式MPC在工程实现中相对技巧性的地方我最初漏掉了绝对幅值约束导致优化器为了满足输出跟踪需求把增量逐步累加后控制量超出执行机构物理上限。4.3 调用quadprog求解并实施第一个增量options optimoptions(quadprog, Display, off, Algorithm, interior-point-convex); [delta_u_opt, fval, exitflag] quadprog(H, f, A_cons, b_ub, [], [], lb, ub, [], options); % 只取第一个增量实施 delta_u_applied delta_u_opt(1:nu); % 更新控制器输出 u_actual u_prev delta_u_applied;这里的u_prev是上一采样时刻实际下发的控制量。实施完第一个Δu后把u_actual保存下来在下一采样时刻作为新的历史值参与状态增广向量拼接。这里有一个容易忽略的小坑quadprog在解半正定问题时会报warning甚至直接退出。增量式MPC的H矩阵理论上是正定的因为R_bar通常取正定对角阵但如果R_bar里某个权重取0H就可能变成半正定此时quadprog内点法照样能解但Matlab新版本会提醒你矩阵奇异。我一般给R_bar对角元素加一个1e-6级别的微小正则项既不影响控制品质又能让求解器稳定运行。4.4 参考轨迹处理前馈还是纯反馈增量式MPC有一个天然特性由于系统内含积分环节即使参考轨迹阶跃变化稳态偏差也会自动归零。但预测时域内如何给参考值直接影响动态响应。我通常分成两档保守档参考轨迹设为恒定目标值不在预测窗口内做规划。优点是鲁棒性强缺点是大阶跃时响应偏慢。激进档参考轨迹按一阶惯性滤波r_filtered(ki) alpha * r_target (1-alpha) * r_filtered(ki-1);alpha越小路径越平滑MPC越容易跟踪alpha越大响应越快但容易触达约束边界。实测下来alpha取0.3-0.5对于大多数过程对象是个不错的中间选择。5. 实测效果与数值病态问题双积分器上的增量式MPC验证5.1 仿真工况设计我在Matlab里用双积分器系统验证上面的公式。采样周期0.1s系统矩阵A [1 0.1; 0 1] B [0; 0.1] C [1 0]目标状态分量x1从0阶跃跟踪到10同时让x2速度在过渡过程中保持相对平稳。参数设置预测时域 Np 30控制时域 Nc 5输出权重 Q 1增量权重 R 0.01Δu约束 ±0.5u绝对幅值约束 ±5完整仿真循环代码结构% 初始化 x [0; 0]; u_prev 0; u_log []; x_log x; ksi [x; u_prev]; for k 1:100 % 构建预测矩阵同上可以预先算好不需要每次重建 F_ksi F * ksi; E r_target - F_ksi; f -Phi * Q_bar * E; % 约束矩阵 u_prev_rep u_prev * ones(Nc, 1); A_cons tril(ones(Nc, 1)); % SISO情况简化 b_ub 5 - u_prev_rep; b_lb -5 - u_prev_rep; % 求解 delta_u quadprog(H, f, A_cons, b_ub, -A_cons, b_lb, -0.5, 0.5, [], options); % 实施 u u_prev delta_u(1); x A * x B * u; % 更新状态 u_prev u; ksi [x; u_prev]; u_log [u_log, u]; x_log [x_log, x]; end实际运行下来最直观的感受是控制量曲线几乎没有抖动。对比位置式MPC同样条件下位置式MPC的u曲线会频繁出现±0.2左右的锯齿形跳动增量式MPC的u曲线则平滑得多只在阶跃开始阶段有一个短暂上升随后非常平稳地趋近稳态值。5.2 数值病态权重矩阵条件数与求解稳定性增量式MPC在实际求解中碰到的最大拦路虎是数值病态。原因是增广后的A_tilde矩阵包含了一个单位块eye(nu)当系统本身有较大极点比如快系统时A_tilde的条件数会变得很大。更糟的是预测矩阵F中A_tilde^i随着i增大会产生巨大的数值量级差异——第1行元素可能是10^0量级第20行元素已经到10^6量级。这直接导致H矩阵的条件数飙升到10^12甚至更多。我处理这类问题有三个实用手段第一无量纲化。这是最根本的解决办法。把输出量、控制量都除以各自的量程让它们在0.1到10之间。比如温控系统输出除以200度阀门指令除以100%开度。无量纲化之后的条件数至少能降两三个数量级。第二缩短控制时域。我实测发现Nc从30降到5H矩阵条件数能降一个数量级以上。背后的道理是Phi矩阵的列数变少Phi*Q_bar*Phi的最小特征值不容易被压到0附近。第三微小正则项。在H矩阵对角线统一加1e-6 * eye(nx_tilde)在很多情况下能救回一个半正定问题。这个方法在quadprog报“Hession matrix is not symmetric positive definite”时特别好用。这三招我在不同项目里反复使用尤其是无量纲化几乎每个模型预测控制落地项目我都会先做归一化处理。它可以避免你在调试时把大量时间浪费在“和求解器作斗争”上让人能更专注地调控制参数。5.3 一个反直觉现象增量权重不能设为零我在调试中试过把增量权重R设成0想看看纯跟踪效果。结果出乎意料系统确实能跟踪但控制量在稳态附近出现难以收敛的高频小抖动而且抖动幅度虽然小幅值执行机构却能明显听到“嗡嗡”声。这个现象的机理是当R0时二次规划的目标函数只惩罚输出偏差而Δu没有任何代价优化器在多个等价解之间任意切换。即使Δu约束存在也只是把切换限制在可行域内不能消除切换本身。这是增量式MPC和位置式MPC的一个微妙差别位置式MPC即使R0输出曲线通常也稳定因为u本身被限制在可行域而增量式中R0相当于对控制动作的频率没有惩罚执行机构长期处于“被优化器反复指挥”的状态。所以我在工程实践中R的初始值从不给0一般先给1e-3量级相对于输出权重归一化后然后逐步减小试探下限直到出现抖动再回调。这个“从大往小试”的方向比“从小往大加”要可靠得多。6. 从双积分器到通用过程模型你的代码能做哪些扩展6.1 输出约束的处理前面实现的代码里没有加入输出约束。实际工程中输出约束几乎不可避免——液位不能超过罐高、温度不能超过材料耐受极限、速度不能超过安全阈值。增量式MPC中加入输出约束的公式是Y F·ξ(k) Φ·ΔU ≤ Y_max转换成关于ΔU的不等式Φ·ΔU ≤ Y_max - F·ξ(k)Matlab左边拼进A_cons矩阵右边拼进b_ub向量。注意输出约束的软硬性质——建议采用软约束形式引入松弛变量否则模型失配时很容易出现优化器无解。我在实际项目里把输出约束分成两段硬约束比如容器不溢出和软约束比如温度尽量不超过某个舒适上限软约束的松弛代价作为额外的线性项加进目标函数。6.2 扰动估计与状态观测器增量式MPC默认模型准确但工程对象谁都不敢保证模型和实际完全一致。标准做法是加状态扰动估计d_hat(k) y_meas(k) - C_tilde·ξ_hat(k)然后在预测时把扰动补偿项加入预测输出Y_pred F·ξ(k) Φ·ΔU d_hat(k)严格来说这属于DMC动态矩阵控制习惯的做法但套在状态空间框架里也完全成立而且实现成本极低——只需要在预测表达式里加上一个常值向量不需要改动H和f的主体结构。这个技巧在很多MPC落地项目里被称为“无模型自适应补偿”我试过对模型偏差30%的对象加了扰动补偿后稳态误差仍然能收敛到1%以内。6.3 多变量系统MIMO的维度变化前面的SISO代码扩展到MIMO只需注意维数匹配nu变成多个输入维度B_tilde变成[(nxnu) × nu]下半块由单位阵扩展而来ΔU向量的长度变成Nc*nu控制量幅值约束的三角矩阵变为分块三角结构每个块是nu×nu的单位阵R_bar需要为每个输入通道单独配置权重而不是用一个标量Matlab里维数一变最容易出错的地方是kron指令的行列安排。我的经验是先画清楚维度表再写代码每拼一个矩阵就用size检查一遍否则很可能浪费一小时在维度不匹配的报错上。% MIMO控制器时域变量维度 % Y: Np*ny × 1 % ΔU: Nc*nu × 1 % Phi: Np*ny × Nc*nu % Q_bar: Np*ny × Np*ny % R_bar: Nc*nu × Nc*nu6.4 从离线仿真到实时部署Matlab仿真跑通之后真正拿到实时系统上跑还有几个代码层面需要注意的事第一预测矩阵可以离线算在线只算f和约束向量。不要每个采样周期都重建F和Phi那样白浪费算力而且代码结构不清爽。把不变的部分抽出来在初始化时算好。第二求解器选项要选对。quadprog的interior-point-convex算法对于中等规模问题比主动集法更快但问题维度如果超过一两百个决策变量建议考虑换用其他求解器比如osqp或嵌入式专门求解器Matlab自带的在实时性上不一定够用。第三时刻记录exitflag和求解时间。我习惯在每个采样周期把这两样东西存进日志变量一旦出现求解失败或者超时事后回放数据时能立刻锁定是哪一步出了问题。7. 从调参到工程验证增量式MPC的实操要点补充7.1 参数整定的先后顺序增量式MPC参数比PID多但整定顺序清晰先固定Np足够大通常取系统上升时间的5到10倍对应的步数再固定Nc为Np的1/5到1/3调Q权重比确定不同输出通道的优先级最后调R增量权重控制执行机构动作剧烈程度我实测中最影响响应品质的是R和Nc的组合R越小、Nc越大响应越快但抖动风险越高反之则响应偏慢但非常稳健。先大步粗调再小步微调整个过程很快就能定位到合适的参数区间。7.2 无模型对象下如何快速判断MPC公式是否正确很多初学者把MPC代码写完后不知道如何判断公式推导和代码实现到底对不对。我给一个特别简单有效的自测方法把MPC当作一维系统测试并与解析最优控制对比。对一阶系统x(k1)a·x(k)b·u(k)在无约束、NpNcN的极限情况下增量式MPC的闭环等效于一个线性状态反馈u(k) u(k-1) - K_fb·[x(k); u(k-1)]把仿真得到的闭环极点和理论极点对比如果一致说明预测矩阵、权重矩阵、求解调用全套逻辑都没问题。这个验证法是我在开发MPC代码后必做的“冒烟测试”一次能排查掉80%的潜在错误。7.3 真实控制对象的采样周期与执行机构响应匹配最后再强调一次采样周期与增量MPC的匹配问题。增量式MPC天然适合采样周期中等偏慢的场景——采样周期太短比如1毫秒而执行机构响应需要几十毫秒优化器算出来的增量往往根本来不及执行浪费算力采样周期太长比如1分钟而系统时间常数只有几秒MPC的预测优势又发挥不出来。我一般按下述经验选采样周期采样周期约为系统主导时间常数的1/10到1/20。比如温度对象时间常数300秒采样周期15秒左右合适伺服系统时间常数0.05秒采样周期2毫秒到5毫秒。在这个区间内增量式MPC的预测价值能充分体现而且不会因为采样太密导致执行机构还没响应完就被下一次命令打断。8. 最后再补一句我自己的体会用输入增量实现状态空间MPC表面上是公式改写本质上是对“控制作用到底是什么”这个问题的重新理解。仿真里用位置式MPC跑得通不代表真实对象能用执行机构对你的控制信号有着天然的“微分响应”特性增量式MPC只是把这个特性变成了控制器公式设计的一部分。代码层面Matlab里的实现核心不外乎三个动作增广状态、构建预测矩阵、拼二次规划问题。每一步都有常见的坑——维度不匹配、权重矩阵病态、约束矩阵符号搞反——但每踩一次坑对MPC的理解就深一层。我自己走了不少弯路才把这套东西跑通所以把完整的公式推导、Matlab代码细节和调试经验整理在这里。如果你正在接触MPC或者被位置式MPC的执行器抖动问题困扰希望这篇东西能帮你省下几个星期的排查时间。后面如果有时间我还会写一篇关于增量式MPC如何与状态观测器结合、以及在Simulink里做硬件在环仿真的文章欢迎继续关注。
