模型预测控制（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)=Ax(k)+Bu(k)\text{（状态方程）}$$
  $$y(k)=Cx(k)+Du(k)\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)=Ax(k)+Bu(k)$$
  $$x(k+2|k)=Ax(k+1|k)+Bu(k+1)=A(Ax(k)+Bu(k))+Bu(k+1)=A^{2}x(k)+ABu(k)+Bu(k+1)$$
  $$\dots$$
  $$x(k+i|k)=A^{i}x(k)+\sum _{j=0}^{i-1}{A}^{i-1-j}Bu(k+j)$$
• 同样可以预测输出 $y(k+i|k)=Cx(k+i|k)+Du(k+i)$ (如果 $D=0$，则为 $Cx(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}\Vert y(k+i|k)-r(k+i){\Vert }_Q^2+\sum _{j=0}^{N_c-1}\Vert u(k+j){\Vert }_R^2+\sum _{j=0}^{N_c-1}\Vert \Delta u(k+j){\Vert }_S^2$$
  其中：
  • $r(k+i)$ 是未来时刻 $k+i$ 的期望输出（设定点或参考轨迹）。
  • $\Vert \text{vector}{\Vert }_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$，求解一个带有约束的优化问题：
  $$\underset{U(k)}{\min }J(U(k))$$
  subject to:
  $$x(k+i+1|k)=Ax(k+i|k)+Bu(k+j)\text{(模型动力学)}$$
  $$y(k+i|k)=Cx(k+i|k)\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$：

1. 获取当前状态 $x(k)$（通过测量或状态估计器，如卡尔曼滤波器）。
2. 设定未来参考轨迹 $r(k+1)\dots r(k+N_p)$。
3. 求解上述优化问题，得到最优控制序列 $U^{\ast }(k)=[u^{\ast }(k{)}^T,u^{\ast }(k+1{)}^T,\dots ,u^{\ast }(k+N_c-1{)}^T{]}^T$。
4. 只将第一个控制输入 $u^{\ast }(k)$ 应用到实际系统中。
5. 等待下一个采样时刻 $k+1$，返回步骤 1。

三、 手动计算示例：揭开 MPC 的面纱

让我们来解一道 MPC 的“数学题”。为了手动计算可行，我们选择一个极其简单的系统和设置。

问题设定：

• 系统模型（一阶离散系统）：
  $$x(k+1)=0.8x(k)+0.5u(k)$$
  $$y(k)=x(k)\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)：

1. 获取当前状态： $x(0)=5$。

2. 预测未来状态 (以 $u(0)$ 和 $u(1)$ 为未知数)：

   • $k+1$ 时刻：
     $$x(1|0)=Ax(0)+Bu(0)=0.8\times 5+0.5u(0)=4+0.5u(0)$$
   • $k+2$ 时刻：
     $$x(2|0)=Ax(1|0)+Bu(1)=0.8(4+0.5u(0))+0.5u(1)=3.2+0.4u(0)+0.5u(1)$$
   • $k+3$ 时刻：
     注意：$N_p=3,N_c=2$。根据常见策略，在控制时域之后 ($j\ge N_c$)，控制输入保持不变。即 $u(2)=u(1)$。
     $$x(3|0)=Ax(2|0)+Bu(2)=0.8(3.2+0.4u(0)+0.5u(1))+0.5u(1)$$
     $$=2.56+0.32u(0)+0.4u(1)+0.5u(1)$$
     $$=2.56+0.32u(0)+0.9u(1)$$

3. 构建目标函数 $J$ (以 $u(0)$ 和 $u(1)$ 为变量)：
   $$J=\sum _{i=1}^{N_p}Q(x(k+i|k)-r{)}^2+\sum _{j=0}^{N_c-1}Ru(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.5u(0){)}^2+(3.2+0.4u(0)+0.5u(1){)}^2+(2.56+0.32u(0)+0.9u(1){)}^2+0.1u(0{)}^2+0.1u(1{)}^2$$

4. 求解优化问题（无约束情况）：
   目标是找到使 $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.5u(0))(0.5)+2(3.2+0.4u(0)+0.5u(1))(0.4)+2(2.56+0.32u(0)+0.9u(1))(0.32)+2(0.1)u(0)$$
     $$=(4+0.5u(0))+(2.56+0.32u(0)+0.4u(1))+(1.6384+0.2048u(0)+0.576u(1))+0.2u(0)$$
     $$=(0.5+0.32+0.2048+0.2)u(0)+(0.4+0.576)u(1)+(4+2.56+1.6384)$$
     $$=1.2248u(0)+0.976u(1)+8.1984$$
   • 计算偏导数 $\frac{\partial J}{\partial u(1)}$：
     $$\frac{\partial J}{\partial u(1)}=2(3.2+0.4u(0)+0.5u(1))(0.5)+2(2.56+0.32u(0)+0.9u(1))(0.9)+2(0.1)u(1)$$
     $$=(3.2+0.4u(0)+0.5u(1))+(4.608+0.576u(0)+1.62u(1))+0.2u(1)$$
     $$=(0.4+0.576)u(0)+(0.5+1.62+0.2)u(1)+(3.2+4.608)$$
     $$=0.976u(0)+2.32u(1)+7.808$$
   • 令偏导数为 0，得到线性方程组：
     $$1.2248u(0)+0.976u(1)=-8.1984(\text{Eq 1})$$
     $$0.976u(0)+2.32u(1)=-7.808(\text{Eq 2})$$
   • 求解这个 $2\times 2$ 线性方程组：
     从 Eq 2 解出 $u(0)$： $u(0)=\frac{-7.808-2.32u(1)}{0.976}$
     代入 Eq 1:
     $$1.2248\left[\frac{-7.808-2.32u(1)}{0.976}\right]+0.976u(1)=-8.1984$$
     $$1.255(-7.808-2.32u(1))+0.976u(1)\approx -8.1984(\text{using }1.2248/0.976\approx 1.255)$$
     $$-9.799-2.912u(1)+0.976u(1)\approx -8.1984$$
     $$-1.936u(1)\approx -8.1984+9.799$$
     $$-1.936u(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$。

5. 考虑约束：

   • 我们计算出的 $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.5u(1){)}^2+(2.56+0.32(-2)+0.9u(1){)}^2+0.1(-2{)}^2+0.1u(1{)}^2$$
     $$J=(4-1{)}^2+(3.2-0.8+0.5u(1){)}^2+(2.56-0.64+0.9u(1){)}^2+0.1\times 4+0.1u(1{)}^2$$
     $$J=3^2+(2.4+0.5u(1){)}^2+(1.92+0.9u(1){)}^2+0.4+0.1u(1{)}^2$$
     $$J=9+(5.76+2.4u(1)+0.25u(1{)}^2)+(3.6864+3.456u(1)+0.81u(1{)}^2)+0.4+0.1u(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.16u(1{)}^2+5.856u(1)+18.8464$$
   • 求导 $\frac{dJ}{du(1)}$ 并令其为 0：
     $$\frac{dJ}{du(1)}=2\times 1.16u(1)+5.856=0$$
     $$2.32u(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^{\ast }(0)=[-2,-2{]}^T$。
   • 注意： 这种逐个钳位的方法在多变量情况下可能不准确。严格的方法是使用支持约束的 QP 求解器。但对于这个简单的例子，这种方法得到了正确结果，因为目标函数单调性明确。

6. 应用第一个控制输入：

   • 根据 MPC 的滚动时域原理，我们只应用计算出的最优序列的第一个元素。
   • 因此，在 $k=0$ 时刻，施加到系统上的控制输入是 $u(0)=-2$。

7. 进入下一时刻 (k=1) - 展望：

   • 系统状态演化：$x(1)=0.8x(0)+0.5u(0)=0.8\times 5+0.5\times (-2)=4-1=3$。
   • 在 $k=1$ 时刻，MPC 控制器将：
     1. 测量（或估计）到新的状态 $x(1)=3$。
     2. 再次进行预测（这次是预测 $x(2|1),x(3|1),x(4|1)$，基于 $x(1)=3$ 和待优化的 $u(1),u(2)$）。
     3. 构建新的目标函数 $J$（可能是相同的形式，但基于新的预测）。
     4. 求解新的优化问题（寻找最优的 $u(1),u(2)$，并考虑约束）。
     5. 得到新的最优序列 $U^{\ast }(1)=[u^{\ast }(1{)}^T,u^{\ast }(2{)}^T{]}^T$。
     6. 应用第一个控制输入 $u^{\ast }(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 不仅仅是一种控制算法，它代表了一种基于模型进行预测、优化和决策的控制思想。它在过程工业、汽车（如自适应巡航、自动驾驶）、机器人、能源管理等领域取得了巨大成功，并且随着计算能力的提升和算法的发展，其应用范围仍在不断扩大。