企业官网建设流程全解析

模型预测控制的教材和论文很多，但刚接触时会发现同一个问题——不同资料里MPC的公式长得都不一样。有的用绝对输入，有的用输入增量，有的把状态量扩了一倍，有的带了干扰估计器。尤其在做Matlab代码实现时，如果不把公式统一，代码根本写不下去。

这篇内容基于一个实际做过的仿真项目，题目就是【MPC】使用输入增量实现了不同的状态空间MPC公式研究（Matlab代码实现）。我会把增量式MPC的原理讲透，把增广状态空间的构造过程一步步拆开，再把Matlab里从离散化到QP求解的完整链路串起来，最后把四种常见公式形态放到一起对比。如果你正在写MPC代码、或者在复现论文里的公式时被“状态空间格式”绕晕，这篇应该能帮你省下不少时间。

1. 从“公式满天飞”到“一个统一视角”

1.1 为什么同一个MPC，资料里的写法各不相同

MPC的核心思想一句话就能说完：在每个采样时刻，基于系统当前状态求解一个有限时域的开环优化问题，只取第一个控制量作用到对象上，下一时刻再重新求解。但同是这个优化问题，不同作者有不同偏好。

第一种是“绝对输入型”。决策变量直接是未来的控制量序列，比如 (u(k), u(k+1), ..., u(k+N_c-1))，优化问题里出现EAST的是输入的绝对值。这种写法最简单，教材里经常用，但有一个先天缺陷：它没有积分作用。如果模型和实际对象之间存在稳态偏差（几乎必然存在），那么即使控制器在设定值处求得的稳态输入与实际需要的输入不一致，系统也会留下稳态误差。

第二种是“增量输入型”。决策变量变成相邻时刻输入的变化量 (\Delta u(k) = u(k) - u(k-1))，本质就是控制量的差分。这样优化器看到的是“下一步动作改多少”，而不是“下一步设成多大”。因为积分器天然包含在差分结构里，稳态误差问题被直接消掉了。这是工业MPC里最常用的形式，原因后面细说。

第三种是“增广状态空间型”。为了让增量形式能套进标准状态空间框架，把状态量扩充一下，比如把([x(k); u(k-1)])或者([\Delta x(k); y(k)])作为新的状态。这个细节看似只是数学技巧，实际操作中影响非常大——它对目标函数矩阵维度、约束矩阵结构、预测方程批量计算方式全都有连锁影响。

第四种是带外部扰动估计或积分器的扩展形式，比如添加阶跃扰动状态、(\delta)干扰模型等。这一般是实际系统有未知扰动或者模型失配时才需要。

这个项目标题里的“使用输入增量实现了不同的状态空间MPC公式”，核心就是用增量的眼光把这些形态统一起来。理解了这一点，再看论文里的任何MPC公式，都能快速还原成代码。

1.2 输入增量到底改变了什么

用生活化类比说，绝对输入型的MPC像一个每步都要重新“摆好”方向盘角度的司机，增量输入型则是只负责“转多少”的司机。后者天然记住上一次动作，控制量变化是平滑的，且很容易对执行机构的变化速率施加约束。

从数学角度看，把决策变量从(u)换成(\Delta u)会带来三个实质性改变。

第一，系统在优化问题的描述中多了一个“累加器”。实际作用于对象的控制输入变成(u(k) = u(k-1) + \Delta u(k))，这就是积分作用。即使求解器给的(\Delta u)是零，上一时刻的输入也会保持，而不是滑向零或某个默认值。

第二，约束表达方式变了。绝对输入型直接约束(u_{\min} \le u(k) \le u_{\max})，挺方便。换成增量型后，未来的绝对输入值是上一时刻输入加上增量序列的累加，所以输入幅值约束必须写成累加和的形式：(u_{\min} - u(k-1) \le \sum_{i=1}^{j} \Delta u(k+i-1) \le u_{\max} - u(k-1))。约束矩阵从简单的对角阵变成下三角累积矩阵。增量本身的约束则仍然简单：(\Delta u_{\min} \le \Delta u(k) \le \Delta u_{\max})。

第三，系统状态方程的形式变了。如果状态量是原始的物理状态(x)，输入是(\Delta u)，那么原来的方程(x(k+1) = A x(k) + B u(k))里面必须把(u(k))展开成(u(k-1)+\Delta u(k))，但(u(k-1))不是原有的状态分量，这就逼迫我们把它添进去。这就是增广状态空间出现的原因。

1.3 这份Matlab代码解决的核心问题

我当时遇到的问题很典型：复现某篇论文里的MPC算法，论文公式里用的状态量是([\Delta x; y])，但在另一篇资料里则要把(u(k- \1))也作为状态。照着第一份公式写好了预测模型，换到另一份时又得重写。后来发现，只要以“输入增量”为主线，重新统一推导一遍，各种形态其实是同一套预测框架下不同坐标系的表达。

Matlab代码里核心做的工作，简单说是三件事：

• 用输入增量作为决策变量，建立增广状态空间模型；
• 批量构建预测方程中的状态传播矩阵与输入响应矩阵；
• 把优化问题装配成标准二次规划，在滚动时域循环中迭代求解。

代码不复杂，但把这三件事理清楚之后，后面所有变体（带扰动模型、带输出权阵、多变量约束）都是往框架里填东西。

2. 增量式状态空间模型与预测方程推导

2.1 离散状态空间模型的基础假设

做MPC仿真前，第一件事就是把连续系统转换到离散域。仿真和实际工程都是数字控制器，全是在离散时间框架下工作。离散化的标准做法是零阶保持器（ZOH）采样，也就是假设控制输入在一个采样周期内保持恒定，这符合绝大多数执行机构的真实行为。

假设原始连续系统是：

[ \dot{x} = A_c x + B_c u,\quad y = C_c x ]

采样周期为(T_s)，采用零阶保持器离散化后得到：

[ x(k+1) = A_d x(k) + B_d u(k) ] [ y(k) = C_d x(k) ]

Matlab里直接一行c2d就解决：sys_d = c2d(sys_c, Ts, 'zoh')。这里有两个容易出问题的细节。一是采样周期的量纲，如果连续模型的单位是秒，而仿真过程中时间向量用毫秒步长，控制周期会混乱，导致预测模型完全错误。二是如果模型里带有时滞环节，c2d不会自动处理纯时延需要考虑输入延迟建模。

在增量式MPC中，这个离散模型的输入(u(k))会被拆成(\Delta u(k))的累加。所以需要一个扩维的状态向量，才能继续用线性矩阵的形式描述整个系统。

2.2 增广状态空间的构造方式

我在代码里实现了两种增广形式，两种本质上等价，但适用场景不同。

第一种形式是“输入记忆增广”。定义增广状态(\tilde{x}(k) = [x(k); u(k-1)])，那么系统方程变成：

[ \begin{bmatrix} x(k+1) \ u(k) \end{bmatrix} = \begin{bmatrix} A_d & B_d \ 0 & I \end{bmatrix} \begin{bmatrix} x(k) \ u(k-1) \end{bmatrix} + \begin{bmatrix} B_d \ I \end{bmatrix} \Delta u(k) ]

输出方程不变。新状态维度是原来的状态维数加上输入维数。这种形式的好处是，绝对输入量(u(k))可以直接从增广状态的第二个分量读出来，约束处理很方便。代价是状态维数变大，计算量相应增加。

第二种形式是“速度型增广”。定义增量状态(\Delta x(k) = x(k) - x(k-1))，同时把输出(y(k))加入状态，形成(\tilde{x}(k) = [\Delta x(k); y(k)])。推导出的增广矩阵是：

[ \begin{bmatrix} \Delta x(k+1) \ y(k+1) \end{bmatrix} = \begin{bmatrix} A_d & 0 \ C_d A_d & I \end{bmatrix} \begin{bmatrix} \Delta x(k) \ y(k) \end{bmatrix} + \begin{bmatrix} B_d \ C_d B_d \end{bmatrix} \Delta u(k) ]

这个形式最大的优势是直接把输出纳入状态，目标函数里对输出加权的部分变得直白，而且状态量本身自带积分效果。这个形式在某些论文里也叫速度形式MPC，特别适合以输出跟踪为主要目标的场景。

我在代码里两种都写了，方便读者对照。不管用哪一形式，关键在于构造新的矩阵(\tilde{A}, \tilde{B}, \tilde{C})之后，后面的推导流程完全一致。也就是说，增广矩阵一旦写完，预测方程可以完全用通用方式生成。

2.3 预测方程批量矩阵的生成逻辑

有了增广状态空间模型后，下一步就是在当前时刻(k)，根据当前状态(\tilde{x}(k))，预测未来(N_p)步的状态轨迹。这一步要求的不是手动写出每个时刻的方程，而是生成批量矩阵。

定义未来控制增量序列：

[ \Delta U = [\Delta u(k), \Delta u(k+1), \ldots, \Delta u(k+N_c-1)]^T ]

其中(N_c)是控制时域。那么在预测时域(N_p)内，状态预测值可以写成：

[ X_{pred} = \Phi \tilde{x}(k) + \Gamma \Delta U ]

矩阵(\Phi)是把系统矩阵的幂次堆叠起来：([\tilde{A}, \tilde{A}^2, \ldots, \tilde{A}^{N_p}]^T)。矩阵(\Gamma)稍微复杂一点，它是脉冲响应系数的下三角堆叠块。简单说，第(j)步预测时，前(j)个控制增量对当前状态都有贡献，贡献系数就是(\tilde{A}^{j-i}\tilde{B})。

要特别提醒的是：实际编写Matlab代码时不要试图用符号推导把(\Phi, \Gamma)写成闭式解，而是用循环或分块赋值来构建。这个做法可靠得多，特别是当系统维度大或(N_p)较大时，符号推导极易出错。

代码实现的核心段是这样的：

function [Phi, Gamma] = buildPredictionMatrices(Ad, Bd, Np, Nc) n = size(Ad, 1); Phi = zeros(Np * n, n); Gamma = zeros(Np * n, Nc * size(Bd, 2)); for i = 1:Np Phi((i-1)*n+1:i*n, :) = Ad^i; for j = 1:min(i, Nc) Gamma((i-1)*n+1:i*n, (j-1)*size(Bd,2)+1:j*size(Bd,2)) = ... Ad^(i-j) * Bd; end end end

这个函数的输出就构成了预测模型的核心。后面不管目标函数怎么改、约束怎么加，只要(\Phi, \Gamma)这两个矩阵是对的，就成功了一大半。我把这个过程也写入博文，就是为了让大家建立这样一个印象：预测模型不是逐时刻递推的循环，而是一次性批量生成。

3. Matlab实现：从数学模型到可运行的QP

3.1 总体代码结构与数据流

整个Matlab程序的流程可以看作一条流水线：

1. 定义系统参数与离散化
2. 构造增广状态空间矩阵
3. 批量构建预测矩阵
4. 生成QP目标函数与约束矩阵
5. 滚动时域迭代求解

代码结构按这个思路组织，每个阶段都独立成函数，便于调试和复用。如果所有逻辑都堆在一个脚本里，仿真跑到一半想检查某个矩阵会被来回打断思路。

我习惯用一个主脚本驱动整个流程，参数集中放在脚本头部，函数模块单独文件保存。例如：

% 主脚本：main_incremental_mpc.m % 参数设置 Ts = 0.1; Np = 10; Nc = 3; % 模型定义 A = [...]; B = [...]; C = [...]; % 离散化 sys_d = c2d(ss(A, B, C, 0), Ts, 'zoh'); Ad = sys_d.A; Bd = sys_d.B; Cd = sys_d.C; % 增广 [A_aug, B_aug, C_aug] = augmentStateSpace(Ad, Bd, Cd, 'input'); % 预测矩阵 [Phi, Gamma] = buildPredictionMatrices(A_aug, B_aug, Np, Nc); % QP装配 [H, f, Aineq, bineq] = assembleQP(Phi, Gamma, Q, R, x0, u_prev, constraints); % 滚动循环 for k = 1:simSteps ... end

这种结构清晰、易扩展，之后想增加软约束或者经济目标函数，只需要在对应函数里做局部修改。

3.2 预测矩阵与目标函数装配的细节

目标函数的常见形式是最小化预测输出与参考轨迹之间的误差，同时惩罚控制增量大小。写成二次型就是：

[ J(\Delta U) = (Y_{ref} - Y_{pred})^T Q_y (Y_{ref} - Y_{pred}) + \Delta U^T R \Delta U ]

把(Y_{pred} = C_a \Phi \tilde{x}(k) + C_a \Gamma \Delta U)代入，整理成标准QP形式：

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

其中：

[ H = \Gamma^T C_a^T Q_y C_a \Gamma + R ] [ f = -(Y_{ref} - C_a \Phi \tilde{x}(k))^T Q_y C_a \Gamma ]

这里有个常见误区：常数项直接扔掉就行，因为不等式约束与求解器只关心二次项和一次项系数。有的资料会把常数项也留下，那是为了输出目标函数值做分析，不影响最优解。

Matlab里装配QP时，我用的求解器是quadprog。这个内置求解器对中小规模MPC完全够用。调用格式是：

[x_opt, fval, exitflag] = quadprog(H, f, Aineq, bineq, [], [], lb, ub, [], options);

如果约束只有简单的输入增量限幅，可以直接用lb和ub；如果有输入幅值约束或者输出约束，就必须走Aineq/bineq全矩阵路线。

3.3 约束矩阵的装配

约束矩阵是MPC代码里最容易写错的地方。因为决策变量(\Delta U)不是实际控制量，所以“输入幅值约束”不能像绝对输入型那样直接赋上下限。

设当前时刻的上一时刻控制输入为(u_{prev} = u(k-1))，未来第(j)步的绝对控制输入为：

[ u(k+j) = u_{prev} + \sum_{i=1}^{j} \Delta u(k+i-1) ]

所以幅值约束写成：

[ u_{min} - u_{prev} \le \sum_{i=1}^{j} \Delta u(k+i-1) \le u_{max} - u_{prev} ]

写成矩阵形式，就是在约束矩阵里放一个下三角全1块。例如控制时域(N_c=3)时，累加矩阵是：

S = [1 0 0; 1 1 0; 1 1 1];

那么幅值约束的Aineq部分是：

[ \begin{bmatrix} S \ -S \end{bmatrix} \Delta U \le \begin{bmatrix} u_{max} - u_{prev} \ -(u_{min} - u_{prev}) \end{bmatrix} ]

增广状态里的状态约束（如果需要对状态变量限幅）则要从预测方程里把相应的输出或状态行抽出来，写成关于(\Delta U)的线性不等式，需要单独构建。

我的经验是：把所有约束矩阵拼装放到一个函数里处理，而不是在滚动循环里反复计算。循环里只更新与当前状态和上次输入有关的那部分，也就是(u_{prev})相关项。这样可以省很多不必要的计算。

3.4 滚动时域控制主循环的实现

滚动时域循环是整个仿真过程的核心。仿真结构一般是：

% 初始化 x = x0; u = u0; % 历史记录 xHist = zeros(nx, simSteps); uHist = zeros(nu, simSteps); for k = 1:simSteps % 获取当前增广状态 x_aug = [x; u_prev]; % 或 [dx; y] % 更新目标函数一次项 f [H, f, Aineq, bineq] = updateQP(...); % 求解QP [dU_opt, ~, flag] = quadprog(H, f, Aineq, bineq, [], [], lb, ub, [], opts); if flag < 0 warning('第%d步QP求解失败', k); end % 取第一个控制增量 du = dU_opt(1:nu); u_new = u_prev + du; % 作用于对象 x_new = Ad * x + Bd * u_new; y_new = Cd * x_new; % 记录并更新 ... end

需要注意，在第一步之前，(u_{prev})通常取为稳态输入或者0。如果系统初始状态离设定点很远，第一步计算出的(\Delta u)可能很大，这可能导致执行机构饱和。实际工程中会配合抗积分饱和和软约束来处理，但在仿真阶段，一个合理的初始稳态点就可以避免这个问题。

这个循环看起来平淡无奇，但很多新手会在这里出错——比如忘了更新(u_{prev})，或者每一步都在重复构建H矩阵没有利用不变的部分。实际上H矩阵只与模型和权重有关，并不依赖当前状态，因此完全可以在循环外提前算好；循环内部只需要重新计算f向量和约束条件右端项。

4. 不同公式形态的横向对比与选型建议

4.1 绝对输入型、增量输入型、增广型的本质区别

表格做一个直观对比，把四种常见形态放在一起看。