企业官网建设流程全解析

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(k+1) = A·x(k) + B·u(k) y(k) = C·x(k) + D·u(k)

标准MPC的做法是，在当前时刻k，预测未来Np步的系统输出，把未来Nc步的控制量u(k+i)当作决策变量，最小化一个形如：

J = Σ ||y(k+i) - r(k+i)||²_Q + Σ ||u(k+i)||²_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)

新系统的状态方程可以推导为：

ξ(k+1) = [ x(k+1) ] = [ A B ] [ x(k) ] + [ B ] Δu(k) [ u(k) ] [ 0 I ] [ u(k-1) ] [ I ]

也就是说：

ξ(k+1) = 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(k+i) - r(k+i)||²_Q + Σ ||Δu(k+i)||²_R

这里的R直接约束的就是执行机构的动作剧烈程度。工程上修改R的物理含义变得极其清晰——R越大，Δu被压得越小，阀门动作越柔和。位置式MPC里R压制的是绝对位置，调参时很难直观理解“R=0.1到底是让阀门开度小还是让阀门动作慢”，而增量式MPC里R=0.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(k+1) = A·Δx(k) + B·Δu(k) Δy(k+1) = C·Δx(k+1) + D·Δu(k+1)

再把y(k)作为增广状态的一部分，得到：

[ Δx(k+1) ] = [ A 0 ][ Δx(k) ] + [ B ] Δu(k) [ y(k+1) ] [ 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 方案对比：什么时候用哪种？

我在自己项目里通常先用方案一快速验证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)*ny+1:i*ny, :) = C_tilde * (A_tilde^i); end % 构建Phi矩阵 for i = 1:Np for j = 1:min(i, Nc) Phi((i-1)*ny+1:i*ny, (j-1)*nu+1: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系统（nu=1）来说，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里某个权重取0，H就可能变成半正定，此时quadprog内点法照样能解，但Matlab新版本会提醒你矩阵奇异。我一般给R_bar对角元素加一个1e-6级别的微小正则项，既不影响控制品质，又能让求解器稳定运行。

4.4 参考轨迹处理：前馈还是纯反馈

增量式MPC有一个天然特性：由于系统内含积分环节，即使参考轨迹阶跃变化，稳态偏差也会自动归零。但预测时域内如何给参考值，直接影响动态响应。我通常分成两档：

• 保守档：参考轨迹设为恒定目标值，不在预测窗口内做规划。优点是鲁棒性强，缺点是大阶跃时响应偏慢。
• 激进档：参考轨迹按一阶惯性滤波：

r_filtered(k+i) = alpha * r_target + (1-alpha) * r_filtered(k+i-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.5
• u绝对幅值约束 ±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降到5，H矩阵条件数能降一个数量级以上。背后的道理是Phi矩阵的列数变少，Phi'*Q_bar*Phi的最小特征值不容易被压到0附近。

第三，微小正则项。在H矩阵对角线统一加1e-6 * eye(nx_tilde)，在很多情况下能救回一个半正定问题。这个方法在quadprog报“Hession matrix is not symmetric positive definite”时特别好用。

这三招我在不同项目里反复使用，尤其是无量纲化，几乎每个模型预测控制落地项目我都会先做归一化处理。它可以避免你在调试时把大量时间浪费在“和求解器作斗争”上，让人能更专注地调控制参数。

5.3 一个反直觉现象：增量权重不能设为零

我在调试中试过把增量权重R设成0，想看看纯跟踪效果。结果出乎意料：系统确实能跟踪，但控制量在稳态附近出现难以收敛的高频小抖动，而且抖动幅度虽然小幅值，执行机构却能明显听到“嗡嗡”声。

这个现象的机理是：当R=0时二次规划的目标函数只惩罚输出偏差，而Δu没有任何代价，优化器在多个等价解之间任意切换。即使Δu约束存在，也只是把切换限制在可行域内，不能消除切换本身。

这是增量式MPC和位置式MPC的一个微妙差别：位置式MPC即使R=0，输出曲线通常也稳定，因为u本身被限制在可行域，而增量式中R=0相当于对控制动作的频率没有惩罚，执行机构长期处于“被优化器反复指挥”的状态。

所以我在工程实践中，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变成[(nx+nu) × 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*nu

6.4 从离线仿真到实时部署

Matlab仿真跑通之后，真正拿到实时系统上跑，还有几个代码层面需要注意的事：

第一，预测矩阵可以离线算，在线只算f和约束向量。不要每个采样周期都重建F和Phi，那样白浪费算力，而且代码结构不清爽。把不变的部分抽出来在初始化时算好。

第二，求解器选项要选对。quadprog的interior-point-convex算法对于中等规模问题比主动集法更快，但问题维度如果超过一两百个决策变量，建议考虑换用其他求解器（比如osqp或嵌入式专门求解器），Matlab自带的在实时性上不一定够用。

第三，时刻记录exitflag和求解时间。我习惯在每个采样周期把这两样东西存进日志变量，一旦出现求解失败或者超时，事后回放数据时能立刻锁定是哪一步出了问题。

7. 从调参到工程验证：增量式MPC的实操要点补充

7.1 参数整定的先后顺序

增量式MPC参数比PID多，但整定顺序清晰：

1. 先固定Np足够大（通常取系统上升时间的5到10倍对应的步数）
2. 再固定Nc为Np的1/5到1/3
3. 调Q权重比，确定不同输出通道的优先级
4. 最后调R增量权重，控制执行机构动作剧烈程度

我实测中最影响响应品质的是R和Nc的组合：R越小、Nc越大，响应越快但抖动风险越高；反之则响应偏慢但非常稳健。先大步粗调，再小步微调，整个过程很快就能定位到合适的参数区间。

7.2 无模型对象下如何快速判断MPC公式是否正确

很多初学者把MPC代码写完后，不知道如何判断公式推导和代码实现到底对不对。我给一个特别简单有效的自测方法：把MPC当作一维系统测试，并与解析最优控制对比。

对一阶系统x(k+1)=a·x(k)+b·u(k)，在无约束、Np=Nc=N的极限情况下，增量式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里做硬件在环仿真的文章，欢迎继续关注。