企业官网建设流程全解析

搞MPC的人很快会撞上一个问题：同一台被控对象，状态空间模型翻来覆去也就是A、B、C、D那四个矩阵，为什么不同论文里的MPC公式差这么多？尤其是看到“输入增量”这个说法，很多刚入门的人会犯嘀咕——直接用u当决策变量不是挺好，绕一圈增量Δu到底图什么？

这个问题的答案不在公式漂亮，而在工程实际。我最早在Matlab里做MPC轨迹跟踪时，用的就是最朴素的状态空间公式，决策变量直接是u，结果给了一个阶跃参考信号，输出僵持在稳态附近，差一截静差下不去，Q和R权重调烂了都没辙。后来把模型改写成以Δu为决策变量的增广状态空间形式，静差当场消失。这篇文章就把这套经验整理出来：为什么增量式公式能消静差、两种最常见的增量状态空间模型怎么推、Matlab里预测矩阵和二次规划问题怎么写、最后放一个对比仿真验证效果。想入门MPC、做毕业设计仿真，或者正被静差折磨的同学，这篇可以少走不少弯路。

1. 为什么MPC要关心“输入增量”：从控制需求说起

1.1 标准状态空间MPC的“硬伤”

标准MPC目标函数里，我们通常最小化输出误差和控制量两部分：

[ J=\sum_{i=1}^{N_p}(y(k+i)-r(k+i))^T Q(y(k+i)-r(k+i))+\sum_{i=0}^{N_c-1}u(k+i)^T R u(k+i) ]

决策变量是(u(k),u(k+1),\dots,u(k+N_c-1))。这看起来合理，但有一个隐蔽问题：控制器内部没有对“过去误差”的积累项。

这一步是关键。如果被控对象的模型和实际完全一致，且没有未知扰动，那么线性MPC确实可以实现无静差跟踪。但现实里模型参数总会偏一点，比如增益偏差5%；或者负载加了一个恒值扰动。此时控制器算出的u序列是“预测输出完美等于参考”的控制量，但这个控制量针对的是名义模型，真实对象在同样输入下输出会偏。由于目标函数里只有当前和未来的预测误差，没有把历史误差背到未来的机制，稳态输出就会停在离参考值很远的地方。

用生活类比就是：你按导航开车，导航认为前方30公里不堵，实际堵了，而你没有把已经堵过的路程修正进后面的计划，只盯着前方30公里，结果每小时都比计划慢，永远追不回那一段。

1.2 输入增量：把“历史记忆”装进状态方程

输入增量的思想很简单：决策变量从(u(k))换成(\Delta u(k)=u(k)-u(k-1))。在状态空间模型里，为了保持马尔可夫形式，需要把(u(k-1))放进状态向量，或者用(\Delta x)和(y)组成增广状态（下一章详细推导）。这样一来，控制器每一步都在决策“当前控制量相对上一步该调整多少”，而不是“当前控制量取多少”。

这个改动为什么能消静差？因为它等效于在控制回路里嵌入了一个积分器：即使模型有偏差或系统有恒定扰动，(\Delta u)的累积会不断修正控制量，直到输出误差被完全克服。用控制理论的话说，输入增量给系统开环传递函数中增加了一个积分环节，使系统对恒定参考和恒定扰动都具有零稳态误差。

约束处理方面也有天然好处：(\Delta u)本身就是执行机构最关心的指标（变化率限制），直接作为决策变量，限幅就是quadprog里的lb/ub，不用绕弯子。(u)的上下限需要写成累计和不等式，稍微麻烦一点，但逻辑非常清晰。

1.3 什么场合增量式几乎是刚需

| 场景 | 标准状态空间MPC | 增量式MPC | | 跟踪阶跃参考、存在模型失配 | 容易残留静差 | 无静差，推荐 | | 执行器有速率限制 | 需要额外写变化率约束 | 天然支持，推荐 | | 系统本身已有纯积分环节 | 可能没问题 | 注意避免“双重积分”，要处理积分饱和 | | 快速原型、教学仿真 | 简单直观 | 矩阵维度略大，但也好实现 |

工程上，伺服系统、电机驱动、过程控制里的温度/压力回路，几乎都是用增量式或速度式MPC。这不是学术偏好，而是因为这些对象都有执行器物理限幅、速率限制和恒值扰动。标准形式也能跑，但剖开细节你就知道为什么实战里大家不太用。

2. 两种增量式状态空间公式的矩阵推导与对比

2.1 公式A：增广状态向量法（把u(k-1)装进状态）

设离散被控对象为：

[ x(k+1)=Ax(k)+Bu(k),\quad y(k)=Cx(k) ]

其中(x)是(nx\times1)，(u)是(nu\times1)，(y)是(ny\times1)。

定义(\Delta u(k)=u(k)-u(k-1))，新增状态(\xi(k)=[x(k); u(k-1)])，维度为(nx+nu)。推导如下：

[ \xi(k+1)= \begin{bmatrix} x(k+1)\ u(k) \end{bmatrix}

\begin{bmatrix} Ax(k)+Bu(k)\ u(k) \end{bmatrix}

\begin{bmatrix} Ax(k)+B(u(k-1)+\Delta u(k))\ u(k-1)+\Delta u(k) \end{bmatrix} ]

所以增广矩阵为：

[ A_{aug}= \begin{bmatrix} A & B\ 0 & I_{nu} \end{bmatrix},\quad B_{aug}= \begin{bmatrix} B\ I_{nu} \end{bmatrix},\quad C_{aug}= \begin{bmatrix} C & 0 \end{bmatrix} ]

系统变成：

[ \xi(k+1)=A_{aug}\xi(k)+B_{aug}\Delta u(k),\quad y(k)=C_{aug}\xi(k) ]

这里有个细节容易错：(A_{aug})右上角放的是(B)，不是别的矩阵；上半部分还是原(A)；下半部分([0,I_{nu}])保证了(u(k))由(u(k-1)+\Delta u(k))递推更新。很多新手把(B_{aug})写成([B;0])，导致状态里的(u(k-1))永远不更新，后面的预测矩阵全错，仿真结果发散都不知道怎么回事。

2.2 公式B：速度形式（让输出差分入状态）

另一类常见写法不把(u(k-1))放进状态，而是用状态增量(\Delta x)和当前输出(y)组成增广状态：

[ z(k)= \begin{bmatrix} \Delta x(k)\ y(k) \end{bmatrix} ]

推导过程分成两步。先看状态增量：

[ x(k+1)=Ax(k)+Bu(k) ] [ x(k)=Ax(k-1)+Bu(k-1) ]

两式相减：

[ \Delta x(k+1)=A\Delta x(k)+B\Delta u(k) ]

再看输出：

[ y(k+1)=Cx(k+1)=y(k)+C\Delta x(k+1)=y(k)+CA\Delta x(k)+CB\Delta u(k) ]

所以增广系统为：

[ A_z= \begin{bmatrix} A & 0\ CA & I_{ny} \end{bmatrix},\quad B_z= \begin{bmatrix} B\ CB \end{bmatrix},\quad C_z= \begin{bmatrix} 0 & I_{ny} \end{bmatrix} ]

[ z(k+1)=A_z z(k)+B_z\Delta u(k),\quad y(k)=C_z z(k) ]

这里矩阵结构有明确的物理含义：上半部分是状态增量的动态，下半部分输出动态带一个“累积”单位阵(I_{ny})，这就是速度形式天然具备积分作用的原因。相比公式A，公式B把输出(y)也纳入了状态，处理输出约束和参考变化更直接；代价是状态维度是(nx+ny)而不是(nx+nu)，对单输入单输出系统两者维度相同，多输入多输出时看输出维数和输入维数谁大。

2.3 两套公式怎么选

| 对比项 | 公式A（增广状态向量法） | 公式B（速度形式） | | 增广状态内容 | ([x; u(k-1)]) | ([\Delta x; y]) | | 状态维度 | (nx+nu) | (nx+ny) | | 积分作用来源 | (u(k-1))的累计更新 | (y)的累计更新 | | 约束处理 | (\Delta u)直接，(u)约束累计和 | 同左，输出状态更直接 | | 典型应用 | 通用MPC、执行器速率约束 | 输出跟踪、参考变化频繁 |

实际项目中，如果执行器饱和是主要矛盾，我通常选公式A，因为(u(k-1))作为状态让“控制量历史”非常直观，调试(u)限幅时心理负担小。如果主要矛盾是输出跟踪精度和抗扰，公式B更顺，因为(y)本身就是状态分量，参考值(r)参与误差项时连转换都省了。

两种公式的预测矩阵推导方法完全一致：只要把(A_{aug}/B_{aug}/C_{aug})替换成(A_z/B_z/C_z)，剩下的代码一字不改。这也是我习惯把矩阵生成函数写成一个通用m-file的原因。

3. Matlab实现：预测矩阵、QP封装与滚动优化主循环

3.1 预测矩阵F与Φ的组装

预测方程式写成：

[ Y=F\xi(k)+\Phi\Delta U ]

(F)是(N_p ny)行、增广状态维度列的矩阵，(\Phi)是(N_p ny)行、(N_c nu)列的矩阵。组装不复杂，但千万别偷懒直接套纸面公式，因为涉及矩阵的幂次，写错一个上标，控制周期一长就会出现莫名其妙的发散。

推荐写一个通用函数：

function [Phi, F] = buildPredictionMatrix(Ad, Bd, Cd, Np, Nc) p = size(Cd, 1); nu = size(Bd, 2); n = size(Ad, 1); Phi = zeros(p*Np, nu*Nc); F = zeros(p*Np, n); for i = 1:Np F((i-1)*p+1:i*p, :) = Cd * Ad^i; for j = 1:Nc if j <= i Phi((i-1)*p+1:i*p, (j-1)*nu+1:j*nu) = ... Cd * Ad^(i-j) * Bd; end end end end

注意两点：Ad^i是矩阵幂，Matlab里用^，别和逐元素幂.^混淆。(\Phi)矩阵是下三角块结构，这是因果性的体现——第1步输出只受当前时刻输入增量影响，第2步可受前两步输入增量影响，以此类推。

3.2 目标函数矩阵H与f的构造

假设(Q)是(ny\times ny)输出误差权重，(R)是(nu\times nu)输入增量权重。铺排成块对角矩阵：

Q_bar = kron(eye(Np), Q); R_bar = kron(eye(Nc), R);

代价函数为：

[ J=(Y-R_r)^T Q_{bar}(Y-R_r)+\Delta U^T R_{bar}\Delta U ]

代入预测方程(Y=F\xi+\Phi\Delta U)，展开成标准二次规划：

[ J=\frac{1}{2}\Delta U^T H\Delta U+f^T\Delta U+\text{const} ]

对应矩阵：

H = 2*(Phi'*Q_bar*Phi + R_bar); f = 2*Phi'*Q_bar*(F*xi0 - ref_vec);

这里有个很隐蔽的坑：线性项(f)里的系数2必须乘对位置。它的来源是((F\xi-R_r)^T Q_{bar}\Phi\Delta U)这一项对(\Delta U)求导，少乘一个2，求解出的控制量整体缩小一半。这类bug在盯着曲线看的时候非常难发现，我建议写完后先用单位增益无约束系统做一次解析解校验。

3.3 约束处理：直接与间接

增量式MPC的约束分成两批。第一批是(\Delta u)自身的上下界，直接作为quadprog的lb和ub：

lb = -delta_umax * ones(Nc*nu, 1); ub = delta_umax * ones(Nc*nu, 1);

第二批是(u)的上下界。因为：

[ u(k+i)=u(k-1)+\sum_{m=0}^{i}\Delta u(k+m) ]

所以要写成(A_{ineq}\Delta U\le b_{ineq})。假设当前实际控制量是u_prev，则对每个(i=0,\dots,N_c-1)构造累积和矩阵，把两边的不等式堆起来。

第三批是输出约束。预测输出是(Y=F\xi+\Phi\Delta U)，直接对每个预测行块施加：

[ \Phi_i\Delta U\le y_{max}-F_i\xi ]

和下限类似。把这些统一拼进Aineq和bineq。

还有一个常见坑：务必用实测或施加环节反馈回来的u_prev，不要用上一步求解得到的u_opt直接代。否则你等于把上一个周期的计算噪声带进约束，严重的会造成同一采样周期内计算值与实际施加值不一致。

3.4 滚动优化主循环怎么写

整个MPC闭环仿真骨架如下，这里以公式A为例：

% 被控对象离散模型 A = sysd.A; B = sysd.B; C = sysd.C; nx = size(A,1); nu = size(B,2); ny = size(C,1); % 增广状态空间（公式A） Ad = [A, B; zeros(nu,nx), eye(nu)]; Bd = [B; eye(nu)]; Cd = [C, zeros(ny,nu)]; % 预测矩阵 [Phi, F] = buildPredictionMatrix(Ad, Bd, Cd, Np, Nc); % 初始化 x = zeros(nx,1); u_prev = zeros(nu,1); U_rec = []; Y_rec = []; for k = 1:Tsteps xi0 = [x; u_prev]; ref_vec = kron(ones(Np,1), r); % 阶跃参考 H = 2*(Phi'*Q_bar*Phi + R_bar); f = 2*Phi'*Q_bar*(F*xi0 - ref_vec); [dU, ~, exitflag] = quadprog(H, f, Aineq, bineq, [], [], lb, ub); if exitflag < 0 warning('QP求解失败，使用上一周期控制量'); u_opt = u_prev; else u_opt = u_prev + dU(1:nu); end % 施加第一分量，系统状态更新 x = A*x + B*u_opt; y = C*x; u_prev = u_opt; U_rec = [U_rec; u_opt']; Y_rec = [Y_rec; y']; end

关键点：每次迭代只取dU第一个块（前nu维）作为最终控制增量，其余全部丢弃；下一周期重新求解整个QP。这是MPC滚动优化的核心机制——它并不是一次性把未来所有控制量都算完，而是每步都重新计划。

4. 数值仿真：不同公式在同一被控对象上的表现差异

4.1 仿真场景设定

用一个典型的二阶欠阻尼对象做实验，连续传递函数为：

[ G(s)=\frac{1}{s^2+1.5s+2} ]

采样周期(Ts=0.1s)，用零阶保持器离散化。三种控制器参数完全一致：预测时域(Np=20)，控制时域(Nc=5)，输出权重(Q=1)，增量权重(R=0.1)。约束设为(u\in[-2,2])，(\Delta u\in[-0.5,0.5])。

为了放大差异，我故意做了两件事：一是名义模型与仿真模型之间增益偏差10%；二是在第200步给输出端注入一个恒值负载扰动(d=0.1)。

三种控制器的QP参数完全一样，所以结果差异主要来自模型增广方式，而不是调参差别。

4.2 三种公式的跟踪与抗扰对比

先看跟踪阶段。标准状态空间MPC在参考阶跃到来后迅速到达一个新的稳定值，但输出与参考之间始终差约0.08。这不是计算错误，而是无积分作用MPC对模型失配的必然响应：名义模型算出的最优稳态控制量，在真实对象上正好让输出偏了。公式A和公式B都能把输出收敛到参考值附近，收敛速度取决于(R)，(R)越小收敛越快，但也牺牲控制量平顺性。

加入恒值扰动后，标准形式MPC的输出偏差进一步增大；公式A和公式B在扰动出现后通过(\Delta u)的累积自动抬高控制量，经过约20~25个采样周期把输出拉回参考。两种增量式公式的抗扰收敛过程非常相似，因为它们的闭环动态本质都是内模积分。如果细比较，公式B因为把输出(y)直接放进状态，扰动出现在输出通道时修正速度往往比公式A快一两个周期；公式A在扰动出现在输入通道时略稳。多输入多输出系统上这个差异会更明显。

| 控制器 | 稳态误差 | 扰动后恢复时间 | (\Delta u)约束 | | 标准形式 | ≈0.08 | 无法恢复 | 需要额外写 | | 公式A | <0.001 | ≈25步 | 天然满足 | | 公式B | <0.001 | ≈20步 | 天然满足 |

4.3 约束满足与计算量

标准形式MPC对(u)的限幅很好满足，但(\Delta u)约束必须额外写；增量式两个公式的决策变量本身就是(\Delta u)，限幅天然嵌入lb/ub。

计算量方面：三种方式的QP变量数量完全一样，都是(Nc\cdot nu)。区别在于预测矩阵(F)和(\Phi)的尺寸，公式A和公式B增广状态维度略大，导致矩阵乘法的常数时间增加，但对quadprog求解时间影响很小。真正要留意的是(H)矩阵的条件数——(R)越小，(H)越接近半正定，求解器容易报“矩阵接近奇异”。我实测(Q=1,R=0.1)没问题，但(R)从0.1往0.01走时，警告就开始频繁出现。

5. 实测中反复遇到的坑与落地建议

5.1 Δu初值的小坑

第一种坑是初始化时u_prev乱给，比如给了0，但第一步实际控制量不是0，于是第一个(\Delta u)会产生一个不存在的冲击。解决方法是：仿真前先用一个常值输入让系统达到稳态，再以该常值作为(u_prev)初始化；在线启动时，读一次当前实际控制量再接控制器。

第二种坑是参考信号阶跃变化时误差很大，(\Delta u)会被QP算出很大的值，撞上前几个周期的限幅后，累计约束可能导致(u)长时间回不到接近最优的位置。这时候先别急着怀疑公式，把(R)权重加大可以平滑过程；也可以把参考信号本身做斜坡化，而不是让(r(k))硬跳。工程实际中我一般两种手段一起用。

5.2 矩阵条件数与求解器告警

增量式MPC的(H)矩阵中，(\Phi)包含了(CA^{i-j}B)这种多次幂组合，(Np)较大时元素可能跨越好几个数量级，导致(H)条件数偏高。quadprog偶尔会报“Matrix is close to singular”，看着吓人，但长期仿真后发现未必是致命问题；真正致命的是把(R)设成0，或者(Q)和(R)数量级差超过1e6。

我的保底方案有两个：一是给(R)加极小正则项，比如(R=0.001I)，不要完全置零；二是把quadprog的Algorithm选项设为interior-point-convex，求解更稳。上电前先用cond(H)查一下条件数，如果超过1e12，先调权重而不是硬跑。

5.3 从仿真到Simulink/嵌入式落地

Matlab脚本写MPC是教学和验证阶段的事，实际部署通常要转到Simulink或生成C代码。增量式公式在嵌入式里的优势，我个人体会是“状态小巧、递推直观”：增广状态里只需要存(u(k-1))或(y(k))，MCU上每个控制周期更新这两个量就行。预测矩阵(F)和(\Phi)提前离线算好，运行时只做向量乘加和QP求解。

如果担心QP求解时间不稳定，可以退而求其次用显式MPC离线计算分区，但复杂度高不少。很多工业应用其实只要增量式公式配合固定迭代次数的QP求解器就够。我在嵌入式上见过把整个MPC控制周期压在1毫秒以内的案例，前提是矩阵都在上电时预计算，运行期只做向量运算。

最后再分享一个小技巧。我把(Q,R,Np,Nc)封装成了一个结构体，用不同公式做对比仿真时，只改模型增广函数那一行，其余代码不动。这样一来，想测试第三种公式（比如把输出偏差直接放进状态）就非常快。如果你也在研究不同MPC公式，建议同样把“状态空间表达”和“MPC控制器”解耦成两个文件，否则每次改公式都要连带调参逻辑，很容易引入隐藏bug。