企业官网建设流程全解析

最近和几位做过程控制的朋友聊到一个很有意思的问题：模型预测控制（Model Predictive Control, MPC）里，为什么很多人不用u做决策变量，非要绕一圈用输入增量Δu？更让人头大的是，不同教材和论文里给出的状态空间MPC公式长得还不一样，有的用状态扩展，有的用输出误差增广，到底哪个能直接搬到Matlab里跑？我花了两周时间，把一套基于输入增量的状态空间MPC从公式推导到Matlab代码完整跑通，顺便把两种主流公式做了横向对比。这篇把中间的所有细节捋一遍，包括推导过程、矩阵构造、quadprog求解，还有我踩过的几个坑，给正在入坑MPC的同学做个参考。

1. 问题引入：为什么MPC要用输入增量公式

1.1 直接控制输入会遇到什么尴尬

如果你刚接触MPC，最自然的想法是把未来N步的控制输入序列u(k)、u(k+1)...u(k+N-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(k+1) = A x(k) + B u(k) y(k) = C x(k) + D u(k)

其中x维度n×1，u维度m×1，y维度p×1。注意D是直接馈通矩阵，物理系统中大多数情况下D=0，但推导时我先保留它，后面会讨论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(k+1) = 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)

写成增广矩阵形式：

ξ(k+1) = [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)

如果假设D=0，输出方程会干净很多，变成y(k) = [C 0] ξ(k)。这也是大多数增量MPC实现里的默认假设。

这里的增广矩阵A_e = [A B; 0 I]，B_e = [B; I]，C_e = [C 0]（D=0时）。维度从原来的n变成n+m，增加的维度正好是输入通道数。这样做的好处是模型结构清晰，可以处理非恒定参考轨迹、可以方便地加入输入幅值约束；坏处是状态里多了一堆u的历史值，噪声环境下估计这些状态需要额外小心。

2.3 误差增广法：差分加跟踪误差

另一种思路是从系统差分出发，得到一组以Δx和输出误差e为状态的模型。先定义两个差分变量：

Δx(k) = x(k) - x(k-1) Δu(k) = u(k) - u(k-1)

对原始状态方程做差分：

Δx(k+1) = x(k+1)-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是恒定参考目标。如果D=0，那么：

e(k+1) = y(k+1)-r = y(k)-r + [y(k+1)-y(k)] = e(k) + C[x(k+1)-x(k)] = e(k) + C[AΔx(k)+BΔu(k)]

写成增广矩阵：

[Δx(k+1); e(k+1)] = [A 0; CA I] [Δx(k); e(k)] + [B; CB] Δu(k)

输出量就是e(k)，即：

e(k) = [0 I] [Δx(k); e(k)]

注意这个模型中状态维度是n+p（状态差分加输出误差），不再需要存u(k-1)。它的含义非常直观：状态的一部分是"状态的增量"，另一部分是"输出对参考的误差"，目标函数里直接惩罚e，控制器天然就是一个积分型控制器。

这种形式在推导预测方程的时候也比较舒服，因为优化问题要跟踪的参考r不直接出现在输出方程里，误差e本身已经是跟踪目标。不过它也有一个限制：要求参考轨迹恒定，或者至少能预先知道参考差分Δr，否则e(k+1)的表达式里要额外减去Δr(k+1)，模型就变复杂了。

2.4 两种公式等价吗

结论是：在相同工况下，两种公式最终闭环行为可以实现一致，因为它们描述的是同一个物理系统，只是状态坐标不同。状态扩展法把历史输入当状态，通用性强；误差增广法把输出误差当状态，直观性更强，更适合恒定参考跟踪。

实际项目里，我一般优先用状态扩展法，因为遇到变参考、扰动和D不为零的情况，不用改模型结构。只有在写学术报告时，才特意用误差增广法展示无静差特性。

3. 不同MPC公式的具体形式与实现要点

3.1 基于状态扩展式的预测方程

无论采用哪种增广模型，MPC的最后一步都要把未来Np步的预测输出表示成当前状态和未来Δu序列的线性函数。以状态扩展式为例，假设D=0，模型为：

ξ(k+1) = A_e ξ(k) + B_e Δu(k) y(k) = C_e ξ(k)

令Np为预测时域，Nc为控制时域，且Nc≤Np，未来控制增量为：

ΔU = [Δu(k); Δu(k+1); ...; Δu(k+Nc-1)] （维度Nc*m×1）

迭代展开：

y(k+1) = C_e A_e ξ(k) + C_e B_e Δu(k) y(k+2) = C_e A_e^2 ξ(k) + C_e A_e B_e Δu(k) + C_e B_e Δu(k+1) ... y(k+Np) = 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(k+Nc-1)

把输出堆成列向量Y = [y(k+1); y(k+2); ...; y(k+Np)]，则：

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(k+i) = u(k-1) + Σ_{j=0}^{i} Δu(k+j)

把所有Nc步的控制量写出来：

U = I_u u(k-1) + L ΔU

其中I_u是Nc*m×1的全1? 准确说，每个控制通道有一个初始值，L是一个由0和1组成的下三角块矩阵。例如单输入Nc=3时：

U = [u(k); u(k+1); u(k+2)] = [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(k+Np)^T P E(k+Np)，其中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)*p+1: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)*p+1:i*p, (j-1)*m+1:j*m) = Ce * Mi; Mi = Ae * Mi; end end end

这个函数假设可以一次性把Np步预测和Nc步控制都算出来，没有考虑Nc后续预测中Δu保持为0的情况。如果你设了控制时域Nc小于预测时域Np，上面Theta矩阵中第j>Nc的列保持为0，表示Nc步之后增量不再变化。使用前先在命令窗口打印size，确认和预期一致。

4.2 二次规划问题装配

有了Gamma和Theta，目标函数可以变成标准的二次型。首先定义输出权重矩阵Q和增量权重矩阵R。如果p=1，Qbar = 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); % 离散模型已经是离散的，直接采样时间Ts=1 % 状态扩展矩阵 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输出系统，p=2，n=2，则A_inc维度4×4，B_inc 4×2，C_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会多出一块。很多教材默认D=0，导致读者在实际系统遇到直接馈通时，一仿真就发现预测输出和真实输出对不上。

我的建议是：如果原始模型D确实不为零，先不要直接套模型，尽量通过合理建模把D消掉。比如传感器和执行器之间如果有纯延迟或滤波器，一般可以把D置零并调整B/C矩阵。如果无法消除，需要在预测方程的Theta里把D的影响叠加上——具体做法是定义一个新的C_tilde = [C D]，并且要特别小心第一行预测输出，因为它会受到当前Δu的直接影响。这时Theta矩阵的第一行对应C_tilde*Be? 注意由于D直接乘Δu，表达式与Be块会有交叉，实际推导时按原始状态方程迭代，不要直接沿用D=0的Theta公式。

5.4 参考轨迹非恒定

状态扩展式可以直接处理随时间变化的参考轨迹，只要把未来的ref序列准备好，塞进E_ref = Gamma*z0 - ref_seq即可。而误差增广式需要参考差分项Δr，模型会变成：

e(k+1) = e(k) + C A Δx(k) + C B Δu(k) - Δr(k+1)

如果你不想改模型，可以对参考信号做前馈补偿，或者直接把参考序列设为零均值偏差。更省事的做法是：非恒定参考路径一律用状态扩展式，恒定参考用误差增广式，各取所长。

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，这套增量状态空间的矩阵构造逻辑仍然是核心的底层框架，不会白学的。