模型预测控制（Model Predictive Control, MPC）手工推导计算
一、 MPC 的哲学与动机：为何需要“预测”与“优化”？
在进入技术细节之前，我们先思考一个根本问题：控制的本质是什么？是根据当前状态（或误差）来决定下一步行动，以期达到某个目标。经典的 PID 控制器就是这种“反应式”控制的典范：它观察当前的误差、误差的积分（历史累积）和误差的微分（未来趋势的初步预测），然后根据预设规则（P、I、D 参数）计算出控制量。
然而，现实世界往往更复杂：
系统惯性与延迟： 控制作用需要时间才能显现效果。当下的最优决策，若不考虑未来的影响，可能导致未来的状态恶化。
约束条件： 执行器（如阀门、电机）有其物理极限（最大/最小开度、最大/最小功率），系统状态（如温度、压力）或输出（如产品浓度）也常常有安全或工艺要求的边界。PID 等传统控制器难以直接、系统性地处理这些约束。
多变量耦合： 许多实际系统是多输入多输出（MIMO）系统，改变一个输入可能影响多个输出，控制回路之间相互影响。解耦控制虽是一种方法，但往往效果有限或设计复杂。
优化目标： 我们往往追求的不仅仅是“稳定”在设定点，还可能希望过程尽可能快、能耗尽可能低、产品质量尽可能高等，这涉及多目标的优化问题。
MPC 的诞生，正是为了应对这些挑战。其核心思想，蕴含着一种“深思熟虑、着眼未来、滚动优化”的哲学：
深思熟虑（基于模型）： 不再是简单的“误差反应”，而是利用对被控对象行为的数学描述（即“模型”），来理解控制输入将如何影响系统未来的状态。模型是 MPC 的基石。
着眼未来（预测时域）： 不仅仅看当前，而是向前“看”一段有限的时间（称为预测时域 \(N_p\) 或 \(P\)）。预测在这段时间内，系统在不同的控制策略下会如何演变。
滚动优化（控制时域与优化）： 在预测时域内，寻找一个最优的控制输入序列（通常只确定未来一小段时间，称为控制时域 \(N_c\) 或 \(M\)，且 \(N_c \le N_p\)），这个序列能使得系统未来的行为（根据模型预测）最符合我们的期望（由一个成本函数或目标函数 \(J\) 来量化），同时满足所有已知的约束条件。
滚动实施（Receding Horizon）： 关键的一步！虽然我们计算出了一段未来的最优控制序列 \(u(k), u(k+1), \dots, u(k+N_c-1)\)，但我们只将这个序列的第一个控制动作 \(u(k)\) 应用于实际系统。然后，在下一个采样时刻 \(k+1\)，我们重新测量（或估计）系统的当前状态，再次进行预测和优化，得到一个新的最优控制序列，并再次只实施第一个动作 \(u(k+1)\)。如此“滚动”向前，不断地用最新的信息修正未来的计划。
这个“滚动优化”策略是 MPC 的精髓所在。它使得 MPC 具备了类似人类决策的特点：基于当前的最佳理解（模型）和对未来的展望（预测），制定一个最优计划（优化），但只执行计划的第一步，然后根据实际情况的变化（新的测量值），不断地重新评估和调整计划。 这使其能够有效应对扰动和模型不确定性。
二、 MPC 的数学构建块：像搭建精密仪器
现在，我们来拆解 MPC 的核心数学部件，理解它们如何协同工作。
系统模型（The Crystal Ball）：
作用： 预测系统未来的状态和输出。
形式： 对于线性系统，常用离散时间状态空间模型：
\[x(k+1) = A x(k) + B u(k) \quad \text{（状态方程）}
\]
\[y(k) = C x(k) + D u(k) \quad \text{（输出方程）}
\]
其中：
\(k\) 是离散时间步长。
\(x(k)\) 是系统在时刻 \(k\) 的状态向量（包含描述系统内部状况的关键变量）。
\(u(k)\) 是系统在时刻 \(k\) 的控制输入向量。
\(y(k)\) 是系统在时刻 \(k\) 的输出向量（我们能测量或关心的变量）。
\(A, B, C, D\) 是描述系统动态特性的矩阵。
关键： 模型的准确性直接影响 MPC 的性能。模型辨识或机理建模是 MPC 应用的第一步。非线性系统则需要非线性模型（NMPC），计算更复杂。
预测（Peering into the Future）：
基于当前时刻 \(k\) 的状态 \(x(k)\) 和未来的一系列控制输入 \(u(k), u(k+1), \dots, u(k+N_p-1)\)，我们可以使用模型迭代预测未来的状态和输出：
\[x(k+1|k) = A x(k) + B u(k)
\]
\[x(k+2|k) = A x(k+1|k) + B u(k+1) = A(A x(k) + B u(k)) + B u(k+1) = A^2 x(k) + A B u(k) + B u(k+1)
\]
\[\dots
\]
\[x(k+i|k) = A^i x(k) + \sum_{j=0}^{i-1} A^{i-1-j} B u(k+j)
\]
同样可以预测输出 \(y(k+i|k) = C x(k+i|k) + D u(k+i)\) (如果 \(D=0\)，则为 \(C x(k+i|k)\))。
注意符号 \(x(k+i|k)\) 表示在时刻 \(k\) 预测的时刻 \(k+i\) 的状态。
控制时域 (\(N_c\))： 我们通常只优化未来 \(N_c\) 个控制输入 \(u(k), \dots, u(k+N_c-1)\)。对于 \(N_c\) 之后的控制输入 (\(j \ge N_c\))，通常假定它们保持不变 \(u(k+j) = u(k+N_c-1)\) 或为 0，以简化计算。
目标函数/成本函数 \(J\) (Defining "Goodness")：
作用： 量化预测的系统行为有多“好”。MPC 的目标是找到使 \(J\) 最小的控制序列。
常见形式（二次规划 Quadratic Programming, QP）：
\[J = \sum_{i=1}^{N_p} \| y(k+i|k) - r(k+i) \|_Q^2 + \sum_{j=0}^{N_c-1} \| u(k+j) \|_R^2 + \sum_{j=0}^{N_c-1} \| \Delta u(k+j) \|_S^2
\]
其中：
\(r(k+i)\) 是未来时刻 \(k+i\) 的期望输出（设定点或参考轨迹）。
\(\| \text{vector} \|_W^2 = \text{vector}^T W \text{vector}\) 是加权二次范数。
第一项：跟踪误差。惩罚预测输出 \(y\) 偏离参考值 \(r\) 的程度。\(Q\) 是半正定权重矩阵，决定了哪些输出的跟踪精度更重要。
第二项：控制能量/幅度。惩罚控制输入 \(u\) 的大小。\(R\) 是正定权重矩阵，用于限制控制动作的剧烈程度，节省能量。
第三项（可选）：控制增量。\(\Delta u(k+j) = u(k+j) - u(k+j-1)\)。惩罚控制输入的快速变化。\(S\) 是半正定权重矩阵，用于平滑控制动作，避免执行器过度磨损。
权重 \(Q, R, S\)： 是 MPC 的关键调优参数。它们反映了不同性能指标之间的权衡（trade-off）。\(Q\) 大则优先保证跟踪精度，\(R\) 大则优先节省控制能量，\(S\) 大则优先保证控制平稳。
约束条件 (Respecting Reality)：
作用： 确保预测和优化的结果在物理或操作允许的范围内。
常见形式：
输入约束：\(u_{\min} \le u(k+j) \le u_{\max}\) (如阀门开度 0-100%)
输入变化率约束：\(\Delta u_{\min} \le \Delta u(k+j) \le \Delta u_{\max}\) (如电机加速度限制)
状态约束：\(x_{\min} \le x(k+i|k) \le x_{\max}\) (如反应器温度限制)
输出约束：\(y_{\min} \le y(k+i|k) \le y_{\max}\) (如产品浓度范围)
关键优势： MPC 能够显式地将这些约束纳入优化问题，这是其相比 PID 等传统方法的核心优势之一。
优化求解器 (The Brain)：
作用： 在每个采样时刻 \(k\)，求解一个带有约束的优化问题：
\[\min_{U(k)} J( U(k) )
\]
subject to:
\[x(k+i+1|k) = A x(k+i|k) + B u(k+j) \quad \text{(模型动力学)}
\]
\[y(k+i|k) = C x(k+i|k) \quad \text{(假设 } D=0 \text{)}
\]
\[\text{Constraints on } u, \Delta u, x, y
\]
其中 \(U(k) = [u(k)^T, u(k+1)^T, \dots, u(k+N_c-1)^T]^T\) 是待优化的控制序列向量。
对于线性模型、二次目标函数、线性约束，这是一个二次规划（QP）问题。有成熟高效的算法（如内点法、活动集法）可以求解。
对于非线性模型（NMPC），这是一个非线性规划（NLP）问题，计算量大得多，求解更困难，可能需要迭代方法（如序列二次规划 SQP）。
计算负担： 优化是 MPC 的主要计算瓶颈，尤其对于快速系统或复杂模型。
滚动时域策略 (The Loop)：
在时刻 \(k\)：
获取当前状态 \(x(k)\)（通过测量或状态估计器，如卡尔曼滤波器）。
设定未来参考轨迹 \(r(k+1) \dots r(k+N_p)\)。
求解上述优化问题，得到最优控制序列 \(U^*(k) = [u^*(k)^T, u^*(k+1)^T, \dots, u^*(k+N_c-1)^T]^T\)。
只将第一个控制输入 \(u^*(k)\) 应用到实际系统中。
等待下一个采样时刻 \(k+1\)，返回步骤 1。
三、 手动计算示例：揭开 MPC 的面纱
让我们来解一道 MPC 的“数学题”。为了手动计算可行，我们选择一个极其简单的系统和设置。
问题设定：
系统模型（一阶离散系统）：
\[x(k+1) = 0.8 x(k) + 0.5 u(k)
\]
\[y(k) = x(k) \quad \text{(状态即可测输出)}
\]
目标： 将系统状态从初始值 \(x(0) = 5\) 调节到设定点 \(r = 0\) (即镇定问题)。
MPC 参数：
预测时域 \(N_p = 3\)
控制时域 \(N_c = 2\)
权重矩阵：\(Q = 1\) (惩罚状态偏离 0)， \(R = 0.1\) (惩罚控制能量) (注意：这里 \(Q\) 和 \(R\) 是标量，因为是 SISO 系统)
约束：控制输入 \(u(k)\) 必须满足 \(-2 \le u(k) \le 2\)。
任务： 计算第一个采样时刻 \(k=0\) 的最优控制输入 \(u(0)\)。
求解步骤 (时刻 k=0)：
获取当前状态： \(x(0) = 5\)。
预测未来状态 (以 \(u(0)\) 和 \(u(1)\) 为未知数)：
\(k+1\) 时刻：
\[x(1|0) = A x(0) + B u(0) = 0.8 \times 5 + 0.5 u(0) = 4 + 0.5 u(0)
\]
\(k+2\) 时刻：
\[x(2|0) = A x(1|0) + B u(1) = 0.8 (4 + 0.5 u(0)) + 0.5 u(1) = 3.2 + 0.4 u(0) + 0.5 u(1)
\]
\(k+3\) 时刻：
注意：\(N_p=3, N_c=2\)。根据常见策略，在控制时域之后 (\(j \ge N_c\))，控制输入保持不变。即 \(u(2) = u(1)\)。
\[x(3|0) = A x(2|0) + B u(2) = 0.8 (3.2 + 0.4 u(0) + 0.5 u(1)) + 0.5 u(1)
\]
\[= 2.56 + 0.32 u(0) + 0.4 u(1) + 0.5 u(1)
\]
\[= 2.56 + 0.32 u(0) + 0.9 u(1)
\]
构建目标函数 \(J\) (以 \(u(0)\) 和 \(u(1)\) 为变量)：
\[J = \sum_{i=1}^{N_p} Q (x(k+i|k) - r)^2 + \sum_{j=0}^{N_c-1} R u(k+j)^2
\]
因为 \(r=0\), \(k=0\), \(N_p=3\), \(N_c=2\), \(Q=1\), \(R=0.1\):
\[J = 1 \cdot x(1|0)^2 + 1 \cdot x(2|0)^2 + 1 \cdot x(3|0)^2 + 0.1 \cdot u(0)^2 + 0.1 \cdot u(1)^2
\]
\[J = (4 + 0.5 u(0))^2 + (3.2 + 0.4 u(0) + 0.5 u(1))^2 + (2.56 + 0.32 u(0) + 0.9 u(1))^2 + 0.1 u(0)^2 + 0.1 u(1)^2
\]
求解优化问题（无约束情况）：
目标是找到使 \(J\) 最小的 \(u(0)\) 和 \(u(1)\)。这是一个关于 \(u(0)\) 和 \(u(1)\) 的二次函数。最小值点出现在偏导数为 0 的地方。
计算偏导数 \(\frac{\partial J}{\partial u(0)}\)：
\[\frac{\partial J}{\partial u(0)} = 2(4 + 0.5 u(0))(0.5) + 2(3.2 + 0.4 u(0) + 0.5 u(1))(0.4) + 2(2.56 + 0.32 u(0) + 0.9 u(1))(0.32) + 2(0.1) u(0)
\]
\[= (4 + 0.5 u(0)) + (2.56 + 0.32 u(0) + 0.4 u(1)) + (1.6384 + 0.2048 u(0) + 0.576 u(1)) + 0.2 u(0)
\]
\[= (0.5 + 0.32 + 0.2048 + 0.2) u(0) + (0.4 + 0.576) u(1) + (4 + 2.56 + 1.6384)
\]
\[= 1.2248 u(0) + 0.976 u(1) + 8.1984
\]
计算偏导数 \(\frac{\partial J}{\partial u(1)}\)：
\[\frac{\partial J}{\partial u(1)} = 2(3.2 + 0.4 u(0) + 0.5 u(1))(0.5) + 2(2.56 + 0.32 u(0) + 0.9 u(1))(0.9) + 2(0.1) u(1)
\]
\[= (3.2 + 0.4 u(0) + 0.5 u(1)) + (4.608 + 0.576 u(0) + 1.62 u(1)) + 0.2 u(1)
\]
\[= (0.4 + 0.576) u(0) + (0.5 + 1.62 + 0.2) u(1) + (3.2 + 4.608)
\]
\[= 0.976 u(0) + 2.32 u(1) + 7.808
\]
令偏导数为 0，得到线性方程组：
\[1.2248 u(0) + 0.976 u(1) = -8.1984 \quad (\text{Eq 1})
\]
\[0.976 u(0) + 2.32 u(1) = -7.808 \quad (\text{Eq 2})
\]
求解这个 \(2 \times 2\) 线性方程组：
从 Eq 2 解出 \(u(0)\)： \(u(0) = \frac{-7.808 - 2.32 u(1)}{0.976}\)
代入 Eq 1:
\[1.2248 \left[ \frac{-7.808 - 2.32 u(1)}{0.976} \right] + 0.976 u(1) = -8.1984
\]
\[1.255 (-7.808 - 2.32 u(1)) + 0.976 u(1) \approx -8.1984 \quad (\text{using } 1.2248/0.976 \approx 1.255)
\]
\[-9.799 - 2.912 u(1) + 0.976 u(1) \approx -8.1984
\]
\[-1.936 u(1) \approx -8.1984 + 9.799
\]
\[-1.936 u(1) \approx 1.6006
\]
\[u(1) \approx -0.8268
\]
将 \(u(1)\) 代回 \(u(0)\) 的表达式：
\[u(0) = \frac{-7.808 - 2.32(-0.8268)}{0.976}
\]
\[u(0) = \frac{-7.808 + 1.918}{0.976}
\]
\[u(0) = \frac{-5.89}{0.976}
\]
\[u(0) \approx -6.035
\]
无约束的最优解是 \(u_{\text{unconstrained}} = [-6.035, -0.8268]^T\)。
考虑约束：
我们计算出的 \(u(0) \approx -6.035\) 违反了约束 \(-2 \le u(k) \le 2\)。\(u(1) \approx -0.8268\) 在约束范围内。
当无约束最优解超出可行域时，最优解通常位于边界上。由于目标函数是二次的（碗状），且可行域是凸的（一个矩形），最优解会是无约束解在可行域上的投影，或者更准确地说，是在违反约束的方向上移动，直到碰到边界。
在这种情况下，\(u(0)\) 违反了下界 -2。因此，我们将 \(u(0)\) 钳位（clamp/saturate）到边界值：\(u(0)_{\text{constrained}} = -2\)。
现在，我们需要重新考虑 \(u(1)\) 的最优值，是在 \(u(0)\) 被固定为 \(-2\) 的前提下。我们应该回到目标函数 \(J\)，将 \(u(0) = -2\) 代入，然后只对 \(u(1)\) 求导并令其为 0。
\[J(u(0)=-2, u(1)) = (4 + 0.5(-2))^2 + (3.2 + 0.4(-2) + 0.5 u(1))^2 + (2.56 + 0.32(-2) + 0.9 u(1))^2 + 0.1(-2)^2 + 0.1 u(1)^2
\]
\[J = (4 - 1)^2 + (3.2 - 0.8 + 0.5 u(1))^2 + (2.56 - 0.64 + 0.9 u(1))^2 + 0.1 \times 4 + 0.1 u(1)^2
\]
\[J = 3^2 + (2.4 + 0.5 u(1))^2 + (1.92 + 0.9 u(1))^2 + 0.4 + 0.1 u(1)^2
\]
\[J = 9 + (5.76 + 2.4 u(1) + 0.25 u(1)^2) + (3.6864 + 3.456 u(1) + 0.81 u(1)^2) + 0.4 + 0.1 u(1)^2
\]
\[J = (0.25 + 0.81 + 0.1) u(1)^2 + (2.4 + 3.456) u(1) + (9 + 5.76 + 3.6864 + 0.4)
\]
\[J = 1.16 u(1)^2 + 5.856 u(1) + 18.8464
\]
求导 \(\frac{dJ}{du(1)}\) 并令其为 0：
\[\frac{dJ}{du(1)} = 2 \times 1.16 u(1) + 5.856 = 0
\]
\[2.32 u(1) = -5.856
\]
\[u(1) = -5.856 / 2.32 \approx -2.524
\]
再次检查约束：计算得到的 \(u(1) \approx -2.524\) 也违反了下界 \(-2\)。
因此，在 \(u(0)\) 已经被固定为 \(-2\) 的情况下，最优的 \(u(1)\) 也应该被钳位到下界：\(u(1)_{\text{constrained}} = -2\)。
所以，在时刻 \(k=0\)，考虑约束后的最优控制序列是 \(U^*(0) = [-2, -2]^T\)。
注意： 这种逐个钳位的方法在多变量情况下可能不准确。严格的方法是使用支持约束的 QP 求解器。但对于这个简单的例子，这种方法得到了正确结果，因为目标函数单调性明确。
应用第一个控制输入：
根据 MPC 的滚动时域原理，我们只应用计算出的最优序列的第一个元素。
因此，在 \(k=0\) 时刻，施加到系统上的控制输入是 \(u(0) = -2\)。
进入下一时刻 (k=1) - 展望：
系统状态演化：\(x(1) = 0.8 x(0) + 0.5 u(0) = 0.8 \times 5 + 0.5 \times (-2) = 4 - 1 = 3\)。
在 \(k=1\) 时刻，MPC 控制器将：
测量（或估计）到新的状态 \(x(1) = 3\)。
再次进行预测（这次是预测 \(x(2|1), x(3|1), x(4|1)\)，基于 \(x(1)=3\) 和待优化的 \(u(1), u(2)\)）。
构建新的目标函数 \(J\)（可能是相同的形式，但基于新的预测）。
求解新的优化问题（寻找最优的 \(u(1), u(2)\)，并考虑约束）。
得到新的最优序列 \(U^*(1) = [u^*(1)^T, u^*(2)^T]^T\)。
应用第一个控制输入 \(u^*(1)\) 到系统中。
这个过程会一直重复下去。
这个手动计算的例子揭示了：
预测是核心： 未来状态是基于模型和假设的未来控制来计算的。
优化是引擎： 通过最小化成本函数来找到“最好”的控制动作组合。
约束是现实： 必须考虑物理和操作限制，这往往会改变无约束的“理想”解。
滚动是策略： 只用第一个动作，持续反馈和重新优化，使 MPC 具有鲁棒性。
计算量： 即使是如此简单的例子，手动计算也相当繁琐。实际 MPC 需要高效的数值优化算法和计算能力。
四、 MPC 的深度思考：超越计算
模型质量的悖论： MPC 的性能高度依赖模型精度。但完美模型不存在。MPC 如何在模型失配下工作？滚动时域的反馈机制是关键，它不断用实际测量来校正预测的起点，一定程度上补偿了模型误差。鲁棒 MPC 则更进一步，明确地在优化中考虑模型不确定性。
计算时延的影响： 优化计算需要时间。在快速动态系统中，当计算结果出来时，系统状态可能已经显著变化，导致计算出的“最优”控制是基于过时信息的。实时性是 MPC 应用的一个重要挑战，催生了显式 MPC（预先离线计算好不同状态区域的最优控制律）等技术。
调参的艺术与科学： \(N_p, N_c, Q, R, S\) 的选择对性能至关重要，但往往缺乏完全系统化的方法，需要经验和反复试验。这更像是一门“艺术”。选择不当可能导致系统性能不佳、振荡甚至不稳定。
\(N_p\) 需要足够长以捕捉系统的主要动态和约束影响，但太长会增加计算负担。
\(N_c\) 影响控制的灵活性和计算量，通常 \(N_c \ll N_p\)。
\(Q\) 和 \(R\) 的相对大小决定了响应速度和控制消耗之间的平衡。
稳定性问题： 经典的 MPC（有限时域）并不自动保证闭环稳定性。虽然滚动时域通常表现稳定，但理论保证需要额外的条件，如终端约束集、终端惩罚函数，或使用无限时域预测（理论上）。
与最优控制的关系： MPC 可以看作是一种近似求解复杂约束最优控制问题的方法。它通过有限时域滚动优化，来逼近（或在特定条件下实现）某种最优性。
MPC 的“智能”： MPC 展现的“智能”来源于其前瞻性（预测）和决策能力（优化）。它不是简单的规则映射，而是基于对系统行为的理解和对目标的追求，在约束下动态地做出决策。
五、 总结：MPC 的力量与挑战
模型预测控制（MPC）是一种先进的控制策略，其核心在于利用系统模型进行未来行为预测，并在满足约束条件下在线优化控制输入序列，然后滚动地实施第一个控制动作。
核心优势：
处理约束的系统性方法： 能直接处理输入、状态和输出约束。
多变量系统控制： 自然地适用于 MIMO 系统，处理变量间的耦合。
优化性能： 可以明确地优化经济或性能指标。
前瞻性： 考虑未来行为，对具有延迟或复杂动态的系统有效。
面临挑战：
模型依赖性： 需要相对准确的系统模型。
计算复杂度： 在线优化可能计算量巨大，对硬件和算法要求高。
参数整定： \(N_p, N_c, Q, R, S\) 等参数的选择需要专业知识和调试。
稳定性保证： 需要特定设计才能严格保证闭环稳定性。
MPC 不仅仅是一种控制算法，它代表了一种基于模型进行预测、优化和决策的控制思想。它在过程工业、汽车（如自适应巡航、自动驾驶）、机器人、能源管理等领域取得了巨大成功，并且随着计算能力的提升和算法的发展，其应用范围仍在不断扩大。