企业官网建设流程全解析

MPC（模型预测控制）这几年在工业过程控制、自动驾驶、机器人运动控制里几乎成了标配。但很多入门的朋友跟我一样，一开始接触的都是最朴素的状态空间MPC：直接以控制输入 u(k) 为决策变量，代价函数里对 u 做惩罚，约束也直接加在 u 上。等项目真落地的时候，发现一个问题：模型有误差或者存在外部常值扰动时，静差消不掉，控制器输出在稳态时总差那么一点。后来我认真啃了采用输入增量（Δu）的MPC公式，发现状态空间模型的推导方式还不止一种，不同公式实现出来的效果、代码复杂度、数值稳定性差别都不小。今天这篇就以Matlab代码实现为主线，把“用输入增量构造状态空间MPC”的几种典型公式掰开揉碎，把我的推导思路、仿真代码和踩坑记录一并分享出来。适合已经会基础MPC推导、想搞清楚增量式实现细节的朋友，也适合正在做Matlab+Simulink控制器开发、打算自己手写MPC求解器的工程师参考。

1. 为什么MPC要引入输入增量：从基本原理讲起

1.1 标准状态空间MPC的基本形式

先把标准离散状态空间模型写在前面，后面所有增量式公式都是从这个模型变过去的。设被控对象为线性时不变离散系统：

x(k+1) = A x(k) + B u(k) y(k) = C x(k) + D u(k)

其中 x∈R^n，u∈R^m，y∈R^p。MPC的核心思想是：在当前时刻 k，利用该模型预测未来 N_p 步的输出或者状态，构造一个带约束的优化问题，求出未来 N_c 步的最优控制序列，然后只把第一步控制量送到被控对象，下一时刻滚动优化。

标准MPC通常直接以 u(k), u(k+1), ..., u(k+N_c-1) 为优化变量，目标函数写成：

J = Σ_{i=1}^{N_p} || ŷ(k+i) - r(k+i) ||Q^2 + Σ{i=0}^{N_c-1} || u(k+i) ||_R^2

约束包括输入幅值约束、输入变化率约束、输出约束等。这种做法在理想模型下没问题，仿真效果很漂亮。但是只要被控对象有未建模动态、模型参数失配或者常值扰动，输出就会存在稳态误差。原因是标准MPC没有把积分作用显式放进控制器结构里，而对无自衡对象或者有扰动的系统，没有积分环节就相当于比例控制，静差是必然的。

1.2 输入增量式MPC的核心动机：消除稳态误差和抗扰动

把控制增量 Δu(k)=u(k)-u(k-1) 作为优化变量，本质上是在控制器内部引入了一个离散积分器。可以这样理解：MPC每步解出来的不是绝对位置的控制量，而是控制量的“变化量”，控制器输出的实际控制量通过累加得到：

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

这样一来，即便模型有常值误差，只要输出偏离设定值，优化问题就会持续产生非零的 Δu(k)，累计之后改变 u(k)，直到误差被校正。有了这个离散积分环节，系统的类型数增加，对阶跃扰动和阶跃给定都能实现无静差跟踪。这跟经典控制里在PID中加积分项、在LQR前面加积分器，原理是相通的。只不过MPC通过把 Δu 放到优化变量里，把积分作用自然地嵌入到了预测模型中，不需要额外设计误差积分状态。

另外，把 Δu 作为优化变量还有一个工程上的好处：可以直接对控制增量施加约束。比如执行机构是伺服电机，给定转速变化率不能超过某个值；或者执行器是液压阀，阀位开度变化速度有限制。这时只需要在优化问题里写 u_min ≤ u ≤ u_max、Δu_min ≤ Δu ≤ Δu_max 即可。如果还是以绝对 u 为变量，变化率约束就得写成相邻两个 u 的差值，形式虽然也能写，但在某些控制器平台上（比如Q3P、FORCES等）不如直接用增量更方便。

1.3 输入增量方案下的两种处理思路

同样是用 Δu 作为优化变量，状态空间模型的写法至少有两条路：

• 扩维状态法：把上一时刻的控制输入 u(k-1) 作为新增状态，构成新状态向量 z(k) = [x(k); u(k-1)]，于是原模型转变成以 Δu 为输入的标准状态空间模型。

• 全增量法：状态和输入一起取增量，即 Δx(k)=x(k)-x(k-1)，Δu(k)=u(k)-u(k-1)，再选择合适的新状态（通常是 [Δx(k); y(k)] 或 [Δx(k); e(k)]）构造增广模型。

两条路最终都能得到形如：

z(k+1) = Ã z(k) + B̃ Δu(k) y(k) = C̃ z(k)

的标准状态空间方程，因此可以统一套用常规的MPC预测矩阵推导，只是在状态维度、矩阵构造细节上有差异。我自己的体会是，扩维状态法直观好理解，适合新手；全增量法在输出跟踪和抗扰动设计上更顺手，尤其是对输出方程里直接含输入项 D 的情况，处理起来更规整。下面两部分分别把这两种公式的推导过程完整走一遍，并给出对应的Matlab矩阵装配代码。

2. 状态空间扩维法：把输入增量变成标准状态空间

2.1 扩维建模的原理推导

假设被控对象模型为不带直接传输项 D 的简化形式，这一步是为了减少干扰，实际中 D 非零时后文也会给出处理方法：

x(k+1) = A x(k) + B u(k) y(k) = C x(k)

定义增广状态 z(k) = [x(k); u(k-1)]，维度为 n+m。因为当前时刻 u(k) = u(k-1) + Δu(k)，所以：

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)

写成矩阵形式：

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

即：

z(k+1) = Ã z(k) + B̃ Δu(k)

其中：

à = [ A, B ; 0, I_m ] B̃ = [ B ; I_m ]

输出方程：

y(k) = C x(k) = [ C, 0 ] z(k)

如果原输出方程中有直接传输项，即 y(k)=C x(k)+D u(k)，代入 u(k)=u(k-1)+Δu(k) 得：

y(k) = [ C, D ] z(k) + D Δu(k)

注意此时输出方程多了一项 D Δu(k)。虽然也可以通过增广输出向量把 D Δu(k) 并进状态，但通常MPC预测时输出表达式中会出现“当前时刻输入增量”对应的矩阵，推导预测方程时需要单独处理，略麻烦。这也是很多人推荐先忽略 D 或者通过模型变换消掉 D 的原因。

2.2 扩维后预测模型的矩阵形式

一旦写出标准状态空间模型 z(k+1)=Ã z(k)+B̃ Δu(k)，接下来的预测推导就和普通MPC完全一致了。设预测时域为 N_p，控制时域为 N_c（并令 N_c ≤ N_p），并假设在 N_c 之后控制增量保持为0。定义从 k 时刻起，未来控制增量向量为：

ΔU = [ Δu(k); Δu(k+1); ...; Δu(k+N_c-1) ] ∈ R^{m N_c}

预测状态序列为：

Z = [ z(k+1); z(k+2); ...; z(k+N_p) ]

则：

Z = Φ z(k) + Ψ ΔU

其中 Φ 是块列矩阵，Ψ 是块下三角矩阵，具体形式跟普通MPC一样：

Φ = [ Ã; Ã^2; ...; Ã^{N_p} ]

Ψ =
[ B̃, 0, 0, ...
à B̃, B̃, 0, ...
Ã^2 B̃, Ã B̃, B̃, ...
...
Ã^{N_p-1} B̃, ..., Ã^{N_p-N_c} B̃ ]

然后预测输出 Y = [ ŷ(k+1); ...; ŷ(k+N_p) ] = C_z Z，其中 C_z = blkdiag(C̃, C̃, ..., C̃)。这里 C̃ 根据扩维状态法输出矩阵而定，若原始输出为 C x，则 C̃ = [ C, 0 ]。

在Matlab中，推荐用循环或者 cell 数组构造 Φ 和 Ψ，不要手动展开大矩阵，否则维度一变就错。后面代码部分我会给出实用函数。

2.3 约束处理与代价函数设计

采用输入增量后，优化变量是 ΔU，需要把原约束统一转成关于 ΔU 的不等式。典型约束有：

• 控制增量约束：Δu_min ≤ Δu(k+i) ≤ Δu_max，这直接是 ΔU 的上下界，处理最简单。
• 控制幅值约束：u_min ≤ u(k+i) ≤ u_max。由于 u(k+i) = u(k-1) + cumsum(Δu(k) ... Δu(k+i))，需要把历史累计关系写成矩阵，例如定义 T_u 为下三角块元素全为 I_m 的矩阵，则 U = repmat(u(k-1), N_p,1) + T_u ΔU_new，其中 ΔU_new 扩展到 N_p 维（N_c之后置0也可以见同构扩展）。上下限约束写成： -T_u ΔU ≤ -u_min + repmat(u(k-1),...) T_u ΔU ≤ u_max - repmat(u(k-1),...)
• 输出约束：y_min ≤ ŷ(k+i) ≤ y_max，直接基于 Y = C_z Φ z(k) + C_z Ψ ΔU 写出线性不等式，本质上也是 ΔU 的线性约束。

代价函数中可以包含输出跟踪误差、输入增量惩罚、控制幅值惩罚等。如果以 ΔU 为变量，目标函数标准形式为：

J = (Y - R)^T Q_y (Y - R) + ΔU^T R_du ΔU + U^T R_u U

其中 R 是参考轨迹向量。由于 U 也是 ΔU 的线性函数，最后能转化为一个标准的凸二次规划（QP）问题。解决QP之后，取 ΔU 的第一个分量 Δu(k)，实际控制量 u(k)=u(k-1)+Δu(k) 输出给被控对象。

这里有个容易忽视的细节：代价函数里的输入增量惩罚矩阵 R_du 和输入幅值惩罚矩阵 R_u 并不是一回事。在一些代码实现中，只惩罚增量，不惩罚幅值，这样控制器可能输出一个很大的偏置量；只惩罚幅值、不惩罚增量，又可能出现剧烈抖动的增量信号。通常的做法是两者都保留，但幅值惩罚项要放很小的权重，或者根本不加入目标函数，而只通过约束限制幅值，具体看工程需求。

3. MATLAB代码实现：增量式MPC控制器搭建全流程

3.1 仿真场景和对象模型

为了对比不同公式，我选了一个经典的二阶开环不稳定对象，加上输入输出模型，形式如下（离散化后）：

A = [1.1 0.3; 0 0.8]; B = [0.1; 0.2]; C = [1 0]; D = 0;

采样周期 T_s=0.1s。这个对象开环极点一个在单位圆外（1.1），用普通比例控制很难稳定，很适合拿MPC演示。仿真中设置设定值 r=1，在 t=5s 加入幅值0.1的常值扰动，考察增量式MPC的静差消除能力。

为什么不直接用更简单的稳定对象？因为在稳定对象上，普通MPC和增量式MPC的静差差别可能不明显，尤其模型精确时基本看不出区别，但换到不稳定对象上，模型失配和扰动的效果一下子就被放大了，更适合暴露问题。

3.2 预测模型构建与Horizon设置

用一个增广函数把原始 A,B,C 转换成扩维模型：

function [At, Bt, Ct] = augmentInputModel(A, B, C, m) % 扩维状态法：z = [x; u_prev]，输入为 du n = size(A,1); At = [A, B; zeros(m, n), eye(m)]; Bt = [B; eye(m)]; Ct = [C, zeros(size(C,1), m)]; end

如果原始模型带 D，则增广输出为 Ct = [C, D]，但同时要注意输出预测时额外加一项 Ddu。为了代码整洁，下面的示例先假设 D=0，实际工程如果 D 不为零，要么对模型做预处理，要么在预测函数里单独加 Ddu 项，我会在后面补一段说明。

然后构造预测矩阵 Ψ 和 Φ。一个简单但清晰的循环写法：

function [Phi, Psi] = buildPredictionMatrices(At, Bt, Ct, Np, Nc) n = size(At,1); p = size(Ct,1); m = size(Bt,2); Phi = zeros(p*Np, n); Psi = zeros(p*Np, m*Nc); % 构造 Phi Pk = At; for i = 1:Np Phi((i-1)*p+1 : i*p, :) = Ct * Pk; Pk = At * Pk; end % 构造 Psi H = zeros(n, m*Nc); H(:, 1:m) = Bt; for i = 1:Np row_start = (i-1)*p + 1; row_end = i*p; for j = 1:min(i, Nc) col_start = (j-1)*m + 1; col_end = j*m; if i == 1 && j == 1 Psi(row_start:row_end, col_start:col_end) = Ct*Bt; else % 计算 At^(i-j) * Bt A_pow = At; for q = 1:(i-j) A_pow = At * A_pow; end Psi(row_start:row_end, col_start:col_end) = Ct * A_pow * Bt; end end end end

这个写法里对每项都重新算 At 的幂，效率不高，但胜在直观，控制工程里时域不会太长（Np 一般不超过20），效率可以接受。追求效率的话可以提前把 At^0,...,At^{Np-1} 存进 cell 数组再调用。

3.3 二次规划求解与闭环仿真代码

增量式MPC的核心在每个采样周期求解一个QP。Matlab环境里最省事的是用旧版quadprog接口，新版也可以用optimoptions指定求解器。先搭出 QP 的标准形式：

min 0.5 * x' * H * x + f' * x s.t. A_ineq * x ≤ b_ineq A_eq * x = b_eq lb ≤ x ≤ ub

其中 x 就是 ΔU。代价函数采用：

J = (Y - R)^T Q_y (Y - R) + ΔU^T R_du ΔU

注意这里为了简洁，暂时不加控制幅值惩罚项，只把幅值约束写进不等式里。展开：

Y = Phi * z + Psi * du_vector

所以：

J = (Psidu - (R - Phiz))' * Qy * (Psidu - (R - Phiz)) + du' * Rdu * du = du' * (Psi'QyPsi + Rdu) * du - 2*(R-Phiz)'QyPsidu + const

常数项不影响最优解，故可令：

H_qp = 2 * (Psi' * Qy * Psi + Rdu) f_qp = -2 * Psi' * Qy * (R - Phi * z)

约束方面，假设只有幅值约束 u_min ≤ u ≤ u_max 和增量约束 du_min ≤ du ≤ du_max。注意 du 向量长度 Nc*m，但实际预测时 Np 可能比 Nc 长，增量在 Nc 之后默认0，所以在构建约束矩阵时，要把 Np 内的控制序列用 Nc 个增量表示：

u_k_i = u_prev + [T_u_cum] * du

其中 T_u_cum 是一个 (Npm) x (Ncm) 的矩阵，第 i 块的累加逻辑是：如果 i ≤ Nc，则第 i 块取前 i 个 du 块之和；如果 i > Nc，则取 Nc 个 du 块之和。Matlab里用 kron 和 tril 构造很方便：

T_cum = zeros(Np*m, Nc*m); for i = 1:Np jmax = min(i, Nc); for j = 1:jmax T_cum((i-1)*m+1:i*m, (j-1)*m+1:j*m) = eye(m); end end % 但注意这里每行块累积和还需要乘以单位阵，这个循环更清楚

更紧凑的写法是：

T_cum = zeros(Np*m, Nc*m); B_eye = repmat({eye(m)}, Np, Nc); for i=1:Np for j=1:min(i,Nc) T_cum((i-1)*m+1:i*m,(j-1)*m+1:j*m)=eye(m); end end

不过还是循环直观。然后幅值约束：

U_full = repmat(u_prev, Np,1) + T_cum * du u_min_full ≤ U_full ≤ u_max_full

即：

T_cum * du ≤ u_max_full - repmat(u_prev, Np,1) -T_cum * du ≤ -u_min_full + repmat(u_prev, Np,1)

增量约束直接落在 du 的上下界。把这些凑成 A_ineq, b_ineq, lb, ub。

完整仿真主程序如下（仅核心部分）：

% 参数设置 Np = 10; Nc = 3; Qy = 10 * eye(Np); % 每个输出均有权重，或针对单输出就标量乘单位阵 Rdu = 0.5 * eye(Nc); u_min = -2; u_max = 2; du_min = -0.5; du_max = 0.5; u_prev = 0; x = zeros(2,1); z = [x; u_prev]; Tsim = 30; Nsim = Tsim / Ts; % 预分配存储 y_hist = zeros(Nsim,1); u_hist = zeros(Nsim,1); options = optimoptions('quadprog','Display','off'); for k = 1:Nsim % 测量或者从对象仿真获取当前状态 x % 参考轨迹 r_k = 1; R = repmat(r_k, Np, 1); % 当前增广状态 z = [x; u_prev]; % 预测矩阵 [Phi, Psi] = buildPredictionMatrices(At, Bt, Ct, Np, Nc); % QP 目标 H_qp = 2 * (Psi' * Qy * Psi + Rdu); f_qp = -2 * Psi' * Qy * (R - Phi * z); % 约束 T_cum = zeros(Np, Nc); for i=1:Np for j=1:min(i,Nc) T_cum(i,j) = 1; % 单输入情况下 m=1，直接标量 end end U_full = repmat(u_prev,Np,1) + T_cum * du; % 这里 du 是符号，实际作为变量 A_ineq = [T_cum; -T_cum]; b_ineq = [repmat(u_max,Np,1) - repmat(u_prev,Np,1); ... -repmat(u_min,Np,1) + repmat(u_prev,Np,1)]; lb = repmat(du_min, Nc, 1); ub = repmat(du_max, Nc, 1); % 求解 du_opt = quadprog(H_qp, f_qp, A_ineq, b_ineq, [], [], lb, ub, [], options); du_k = du_opt(1); u_k = u_prev + du_k; % 仿真对象 x_next = A * x + B * u_k; y_k = C * x_next; % 注意输出取 x(k+1) 或者 x(k)+噪声均可，这里简化为下一拍输出 % 或者按 y=Cx 当前状态，按需处理 % 更新 x = x_next; u_prev = u_k; y_hist(k) = y_k; u_hist(k) = u_k; end

上面代码为了示意，做了一些简化（比如输出取的是下一拍状态），但整体闭环逻辑完整。实际部署时建议把预测矩阵的构建移到循环外，因为模型是定常的，Phi 和 Psi 可以提前算好，循环里只更新 z 和参考轨迹，这样仿真速度能快不少。

3.4 结果可视化分析

仿真结束后，画出输出曲线和控制曲线。增量式MPC在 t=5s 常值扰动加入后，输出经过短暂波动能重新回到设定值，这说明积分作用确实生效了。对比普通MPC（以 u 为变量）的代码，两者在初始响应阶段差别不大，但扰动阶段普通MPC的输出会固定在某个偏差处，怎么调 Q、R 都压不掉那点误差，原因就是没有积分环节。

另外我还测试了不同的 Rdu 对系统动态的影响。Rdu 越大，控制增量被约束得越保守，系统响应变慢，但控制量变化更平滑；Rdu 过小时，控制器对误差反应过于激进，可能出现振荡。直观上，Rdu 就类比对控制量的“阻尼系数”，调参时可以先从 Rdu=1 开始，再往两个方向扫。

4. 不同状态空间MPC公式的对比与选型

4.1 增量式输出MPC与增量式状态MPC的差异

全增量法用的是另一种扩维方式，将状态增量 Δ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)

同时：

y(k+1) = C x(k+1) = C A x(k) + C B u(k)

要做到用 Δx 和 y 表示 y(k+1)，需要对 x(k) 做变换。常见技巧是引入输出方程：y(k)=C x(k)，于是 x(k) 可以由当前状态得到。但更标准的做法是定义新状态 z = [Δx(k); y(k)]，并利用以下关系：

y(k+1) = y(k) + C Δx(k+1) = y(k) + C (A Δx(k) + B Δu(k)) = C A Δx(k) + y(k) + C B Δu(k)

不对，上述推导有问题，因为 y(k)=C x(k)，Δx(k+1)=x(k+1)-x(k)，则 y(k+1)=C x(k+1)=C x(k)+C Δx(k+1)=C x(k)+C A Δx(k)+C B Δu(k)。而 C x(k) 不等于 y(k)+C Δx(k)? 实际上 x(k)=x(k-1)+Δx(k)，所以 C x(k)=C x(k-1)+C Δx(k)=y(k-1)+C Δx(k)，这里多出来 y(k-1)。因此标准做法中需要把状态扩成 [Δx(k); y(k)] 的时候，推导利用 x(k)=Δx(k)+x(k-1)，且 Δx(k+1)=A Δx(k)+B Δu(k)。然后：

y(k+1)=C A x(k)+C B u(k) 这个表达式里含绝对状态和绝对输入，不好直接转化。王鹏教材里采用的是对状态方程先取增量，再取输出为 y(k+1)=C A x(k)+C B u(k)，最后巧妙利用 y(k)=C x(k) 以及 x(k)=Δx(k)+x(k-1) 得到：

y(k+1)-y(k)=C A Δx(k)+C B Δu(k) 吗？我们验证：

y(k+1)=C x(k+1)=C(A x(k)+B u(k)) y(k)=C x(k) 因此 y(k+1)-y(k)=C(A-I)x(k)+C B u(k)。包含绝对量，不简洁。实际上需用扩大状态的等式：

取 z1(k)=Δx(k)，z2(k)=y(k) 则 Δx(k+1)=A Δx(k)+B Δu(k) 同时 y(k+1)=y(k)+C Δx(k+1)? 不对，因为 y(k+1)-y(k)=C x(k+1)-C x(k)=C(x(k+1)-x(k))=C Δx(k+1)。 这里注意 C Δx(k+1)=C(x(k+1)-x(k)) = y(k+1)-y(k)。所以 y(k+1)=y(k)+C Δx(k+1)。而 Δx(k+1)=A Δx(k)+B Δu(k)。因此：

y(k+1)=y(k)+C(A Δx(k)+B Δu(k)) = C A Δx(k)+ y(k)+ C B Δu(k)。推导成立。需要利用 y(k)=C x(k) 吗？不需要，只用到 y(k) 作为积分器状态。所以全增量模型为：

[ Δx(k+1) ] [ A 0 ] [ Δx(k) ] [ B ] [ y(k+1) ] = [ C A I ] [ y(k) ] + [ C B ] Δu(k)

输出为 y(k) (其实 y(k+1) 是状态第二分量)。注意这里 y(k) 是“上一时刻输出”的积分，模型内蕴含积分环节，状态维度为 n+p，比扩维状态法的 n+m 可能不同。

这个方法在输出跟踪问题中很直接：新增的状态就是被控输出本身，你可以直接在代价函数里惩罚 y(k+i)-r，而不需要额外从预测输出向量中提取。它本质上把输出积分状态显式建进了模型，而不是靠外部累加 u。这跟输入增量扩维法（状态里放 u_prev）是不同的思想：一个把执行机构位置当作状态，一个把输出误差的积分当作状态。两者都有积分作用，但数值特性、状态可观测性和对应的观测器设计差别很大。

4.2 观测器设计对公式选择的影响

实际系统状态往往不能全测，需要设计状态观测器。这时增量式MPC的状态构成会直接影响观测器设计。

扩维状态法里包含 u(k-1)，这个信号是控制器已知的，完全可以当作已知输入处理，所以观测器只需观测原始 x。若把 z 整体交给观测器，其中 u_prev 已经是精确已知值，等于增加了一个虚拟的“无噪声状态”，观测器设计简便。

全增量法的状态包含 Δx 和 y。其中 y 通常是传感器直接测到的输出，可作为已知量注入观测器；Δx 则需要估计。这时候状态方程里出现了 y(k) 到下一状态 y(k+1) 的传递，观测器矩阵不再是标准形式，需要仔细配置。如果状态维度低、系统阶次小，两者都能用；但换成高阶系统，我建议优先考虑扩维法，工程实现上更容易处理。

4.3 计算量与数值稳定性分析

从矩阵维度看，扩维法增广状态维度是 n+m，全增量法增广状态维度是 n+p。一般 m 和 p 都远小于 n，所以计算量差异不大。但要注意：如果输出数量 p 很大（比如有多路输出），全增量法状态维度会明显增加，QP 的预测矩阵 Ψ 会变胖，求解变慢。而扩维法状态维度取决于输入数量 m，跟输出数无关。因此多输出系统用扩维法更有优势。

数值稳定性方面，扩维法把 u_prev 作为状态，状态矩阵中会出现单位阵块，特征值稳定；全增量法中由于新增输出积分状态，系统矩阵会有位于单位圆上的特征根（离散积分器）。MPC 预测时，这个积分状态可能因为模型失配导致缓慢漂移，尤其是存在测量噪声时，输出积分状态可能累计噪声，造成控制量漂移。解决方法是给积分状态设置抗饱和机制，或者在Q中降低对积分状态初值的敏感度。这些细节教科书里很少提，但工程上非常关键。

5. 调参与踩坑记录：增量式MPC常见问题速查

为了让你少走弯路，我把调参和写代码时撞过的坑整理成速查表：

5.1 控制器初始化和冷启动问题

增量式MPC启动时，u_prev 怎么设很关键。如果对象本来就运行在一个稳态，u_prev 应该用当前实际控制量，而不是默认0。我一开始用0，结果第一步 Δu 需要从0跳变到稳态控制量附近，如果稳恒控制量比较大，直接触发 u 幅值约束，导致初始响应出现奇怪的非线性。正确做法是：启动控制器前先记录一次执行机构反馈值，作为 u_prev 初始值。

5.2 预测时域与控制时域的配合

增量式MPC对Horizon的敏感程度高于普通MPC。Nc 太小时，控制自由度不足，积分作用发挥不出来；Nc 太大，QP求解变慢，且过大的 Nc 可能造成过度调节。我的经验是 Nc 选在系统主导时间常数对应采样步数的10%~20%，Np 选 Nc 的3~5倍。比如采样周期0.1s，系统时间常数约1s，那么 Np 取15~20，Nc 取3~5。这个组合鲁棒性比较好。

5.3 权重矩阵调节技巧

Qy 决定输出跟踪能力，Rdu 决定执行机构活动强度。一个直观思路：先把 Rdu 设为0，看系统最快动态下的控制量是否可接受；然后逐渐增大 Rdu，直到控制量变化幅度降到一个顺眼的水平。实际运行中，若稳态时有轻微噪声，可以在输出误差项基础上加入一个对输出变化率的惩罚，抑制噪声传到控制端，效果比单纯调大 Rdu 更细腻。

5.4 增广状态导致的可观性与数值问题

扩维法状态矩阵中出现了单位阵块，模型本身是可控的，但也可能存在数值病态。比如 A 中某些元素数量级差异很大，计算预测矩阵时会产生大数吃小数。建议在构造预测矩阵前先对模型做离散化尺度归一化，或者用平衡模型（balreal）预处理。全增量法的积分状态在数值上容易累积误差，连续长时间运行时要考虑定期重置积分状态或增加状态估计修正。

5.5 MATLAB实现中的典型Bug

我写代码时遇到过几个经典bug：

• 用 kron(eye(Np), C) 构造输出块对角矩阵时，C 的维度搞错，导致 Phi 维度不匹配。
• 在约束矩阵 T_cum 中把 Np 和 Nc 搞混，累积矩阵成了 Nc 行，程序还能跑，但实际控制序列不对，且难以察觉。
• 把 quadprog 的 H_qp 少乘了2，目标函数等效问题不大但最优解受约束影响时会有偏差。
• 仿真时输出取 x(k) 还是 x(k+1) 不一致，导致延迟一拍，控制器误以为模型有滞后。
• 忘了在每次循环里更新 u_prev，导致控制器输出的不是累计值，而是纯增量，这会让系统变成纯积分控制，极易发散。

针对最后一个坑，我后来会把“更新 u_prev”这行代码放在仿真函数的前面，并用 assert 检查 u 的连续性。

6. 后续扩展思路：非线性系统、自适应MPC等

6.1 从线性MPC到线性时变MPC

增量式状态空间公式很容易扩展到线性时变系统（LTV-MPC），只需要在每个采样时刻更新 A(k), B(k) 并重新构造预测矩阵。比如对非线性系统在工作点线性化，得到 A_lin(k), B_lin(k)，然后用同样的增量公式生成扩维模型。此时预测矩阵 Ψ 的构造依赖各步的时序模型，不能再用不变矩阵，Matlab里可以逐步组装，但计算量明显上升。我建议用 MATLAB function 块或者生成 C 代码后部署到实时仿真器里。

6.2 与无模型方法结合的实现方向

增量式MPC中的 Δu 思想也可以迁移到无模型自适应控制、数据驱动预测控制中。例如不再显式建模 A,B，而是通过历史数据估计预测矩阵，或者直接学习增量模型 f(Δu, Δx)，本质上仍然是利用“增量”作为决策变量来减少模型误差带来的静差。我最近在尝试用神经网络拟合被控对象的增量动态，再把预测控制问题套进去，思路类似，但计算负担更重，适合离线训练加在线快速推理。

最后再分享一个我个人的小习惯：做MPC仿真时，把不同公式的控制器封装成统一的接口，输入 A,B,C,D、时域参数、权重、约束，输出控制量。这样后续换模型、换控制器，只需要替换底层预测模型函数，不用改动主循环。我最初为了比较扩维法和全增量法，写了两套完整循环，结果调参时要改两处，对比起来非常痛苦。后来重构成接口形式，工作量减少了一半，也更容易定位差异到底是公式本身还是代码实现造成的。

这个项目做完，我最大的体会是：MPC公式看似复杂，但核心就是“怎么把积分作用装进预测模型”。输入增量法给出的答案很优雅：把控制增量当成决策变量，状态空间自然扩维。只要吃透了扩维这一步，后面所有推导和代码都是水到渠成。希望这篇分享能帮你跨过这个坎，早日造出自己顺手的状态空间增量MPC控制器。