企业官网建设流程全解析

MPC（模型预测控制）里，优化变量到底是选 u 还是选 Δu，这个问题我在刚入门时被绕了很久。后来在某伺服系统项目里做位置跟踪，被稳态误差和教育了一回，才真正意识到“输入增量”这四个字的含金量。现在回头看，很多讲MPC的资料直接把增广状态空间公式甩出来，却没有把“为什么要把输入变量改成增量”“不同公式之间到底差在哪”讲透。这篇文章就用Matlab代码实现的方式，把我用输入增量构造不同的状态空间MPC公式的过程、推导步骤和踩坑记录整理出来，主要面向正在自己手写MPC、不想只依赖现成工具箱的开发者。

先说清楚范围：我们只讨论离散线性状态空间模型，目标是通过输入增量（Δu）作为优化变量，实现几种不同的MPC公式，包括直接在原状态空间上做增量优化、把输入扩成增广状态、以及输出增量型的velocity form。这三种公式最后都能化成同一个标准QP问题，但预测矩阵、状态初始化和约束处理各有差异，这正好对应标题里说的“不同的状态空间MPC公式”。我会把每一步的Matlab关键代码贴出来，并给出我在仿真中踩过的坑和调参经验。

1. 为什么MPC要用输入增量代替绝对输入

先回到最基础的问题：MPC在每个采样周期都要根据当前测量值，用模型预测未来 N 步，求解一个带约束的优化问题，然后把得到的第一拍控制量送给被控对象。这里的优化变量如果取成绝对输入序列 u_k, u_{k+1}, ..., u_{k+N-1}，就是我们常见的标准MPC形式；取成增量序列 Δu_k, Δu_{k+1}, ..., Δu_{k+N-1}，就变成了增量式MPC。

从控制本质上看，增量式 MPC 相当于在控制回路中自然引入了一个积分环节。原因很简单：实际执行的控制量为 u_k = u_{k-1} + Δu_k，控制器每次决策的是“在前一拍基础上加多少”。即使模型存在稳态误差或者有常值扰动，只要闭环稳定，积分作用会把输出慢慢拉回参考值。如果直接用绝对输入 u 作为优化变量，相当于一个纯比例式的前馈控制方案，遇到模型失配时很容易留下残余偏差。我用一个简单的类比：普通MPC像你开车时只盯着当前位置和目标的距离打方向，增量MPC则额外记住了方向盘上一拍在哪，让你每一拍知道动了多少，不会因为路面变化而越跑越偏。

另外，增量形式还能很自然地处理执行机构的速率约束。伺服电机、液压阀、舵机这些执行器通常都有自己的最大变化率限制，也就是 |Δu| ≤ Δu_max。优化变量本身就是Δu，那这个约束就直接写成变量的上下界，处理起来非常方便。如果优化变量是绝对输入u，你还需要单独再算相邻控制量之差，麻烦不说，还容易把预测方程搞得一团糟。

但增量形式也有代价：状态空间模型必须跟着变。原本的模型是 x_{k+1} = A x_k + B u_k，输入是 u，现在输入变成了 Δu，状态方程不会自动兼容。这时候就得对模型做改造，可以是把上一拍输入 u_{k-1} 塞进状态向量，也可以是用增量状态 Δx 和输出 y 组成一个新的增广系统。这就引出了多种所谓的“状态空间MPC公式”，下面一章就逐个推导。

2. 三种常见状态空间MPC公式推导

为了方便说明，先给定同一个离散线性定常系统：

x(k+1)=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_{k+1}; ...; Δu_{k+N-1}]

实际施加到系统的输入序列：

U = [u_k; u_{k+1}; ...; u_{k+N-1}]

根据 u_{k+j} = u_{k+j-1} + Δu_{k+j}，可以得到：

U = M u_{k-1} + L ΔU

其中 M 是长度为 N 的全1列向量，L 是 N×N 的下三角全1矩阵。比如 N=3 时：

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_{k+1} = A x_k + B u_k = 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

写成矩阵形式：

v_{k+1} = [[A, B], [0, 1]] v_k + [[B], [1]] Δu_k

y_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_x+1，增广矩阵传递函数多了一个积分极点，要求原系统可观测且增广系统可控，否则会出现不可镇定模式，后面仿真时要注意。

2.3 类型C：完全增量状态空间模型（velocity form）

第三种公式更彻底，不只把输入增量，连状态增量 Δx 也纳入新的状态向量，构成所谓 velocity form。定义：

w_k = [Δx_k; y_k]

其中 Δx_k = x_k - x_{k-1}。推导过程稍微绕一点：

Δx_{k+1} = x_{k+1} - 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_{k+1} = C x_{k+1} = C (x_k + Δx_{k+1}) = C x_k + C(A Δx_k + B Δu_k)

因为 y_k = C x_k，所以：

y_{k+1} = y_k + C A Δx_k + C B Δu_k

于是增广状态方程为：

w_{k+1} = [[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 z

s.t. A_eq z = b_eq

A_in z ≤ b_in

lb ≤ 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)*ny+1:j*ny, :) = C * A^j; for i = 1:j Gamma((j-1)*ny+1:j*ny, i) = C * A^(j-i) * B; end end end

注意这里的 F 中用的是 A^j 而不是 A^(j-1)，因为预测的第一行对应 k+1 时刻：y_{k+1} = C(A x_k + B u_k)。Gamma 的第 i 列对应第 i 个输入分量 u_{k+i-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 * L，bias_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或C，u_{k+j} = u_last + Σ_{i=0}^{j} Δu_{k+i}，依然可以用 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)

采样周期 Ts=0.1s，零阶保持器离散化。预测时域 N=20，控制时域 N也用20（先不搞Nc<Np的复杂情况）。权重 Q=10，R=1，输入幅值约束 |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 画出来对比，你会发现一个很有意思的现象：在模型完全准确、没有扰动的情况下，三条响应曲线几乎重合；一旦我在状态更新里注入一个常值扰动 d=0.1，类型A的输出会出现一个固定偏差，而类型B和类型C凭借内部积分作用能把输出重新拉回参考值。这个差别就是“输入增量”的直观体现。

4.3 三种公式的对比结果

我实际仿真时的数据大致如下，但不同对象和权重下会有差别，关键看趋势：

类型C在实际使用中通常是最稳的，因为输出被做成状态后，即使预测模型存在小幅误差，积分作用也会持续修正，不会让输出长期漂移。类型B在数学上等效于类型A加了一个积分状态，但预测矩阵更整齐。类型A虽然代码量也不复杂，但容易在累积矩阵上犯错，而且遇到模型失配时天然吃亏。

5. 常见问题与调参经验

自己在Matlab里手写增量MPC时，我踩过不少坑，这里挑几个最典型的说。

5.1 quadprog报错：Hessian不是正定矩阵

增量MPC的H矩阵原始形式是 G'QbarG + Rbar。如果G的列线性相关，或者Rbar很小，会导致H奇异。最简单的做法是在H上加一个很小的单位阵：

H = H + 1e-8 * eye(N);

这样能保证quadprog有解。但不要加太大的值，否则等于偷偷把增量权重调大，响应会变慢。另外，如果 Q=0，输出误差完全不惩罚，H可能只由Rbar决定，虽然正定但会失去跟踪能力，这类低级错误也要排查。

5.2 类型C起步跳变

类型C的初始状态是 w0 = [x - x_prev; y]。仿真第一拍如果我把 x_prev 设成了0，而 x 是初始状态，那 Δx 就是一个很大的非零向量。这个非零状态经过预测矩阵放大后，会让第一拍增量 Δu 输出一个离谱的值。解决办法是：如果系统初始确实静止在某个工作点，直接令 x_prev = x，令 Δx=0。很多人会犯这个错，因为直觉上“初始状态为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，导致预测输出少了一拍。一个快速检查方法：设 N=1，此时 F 应该是 CA，Gamma 应该是 C*B。如果输出不对，大概率是索引问题。另外，当 Nc<Np 时，控制时域之后的 Δu 要补零，L矩阵的列数小于行数，很多初学者在这里维度报错。我一般建议先让 Nc=Np 跑通，再去优化控制时域。

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 在拼命挣扎但输出没有改善，大概率是预测方向反了或者状态初始化有误。这个检查习惯帮我省下了很多排查时间。