MPC 学习路径

💡 0. 核心理念

为什么很多人学不透 MPC？

常见的错误学习方式：

• 一上来就读 Rawlings 的理论大部头，被稳定性证明劝退
• 只用 Matlab MPC Toolbox 调参，不理解底层 QP 是怎么构造的
• 跳过 LQR 直接学 MPC，不理解代价函数和 Riccati 方程的来龙去脉

四层递进法

1 第一层：基础原理（2-3 周）

在接触 MPC 之前，必须先掌握以下三个基础工具。MPC 本质上是 LQR 的在线滚动版本，没有 LQR 基础就去学 MPC，代价函数的由来永远搞不清楚。

1.1 状态空间模型

MPC 使用状态空间模型进行预测，这是最核心的数学工具。

必须掌握：

• 从微分方程/传递函数到状态空间的转换
• 能控性 (Controllability) 与能观性 (Observability)
• 状态空间模型的物理含义（每个状态变量代表什么）

1.2 LQR 最优控制

LQR 是 MPC 代价函数的直接前身。理解 LQR = 理解 MPC 代价函数的灵魂。

LQR 代价函数：

LQR 解：

其中 \(P\) 满足代数 Riccati 方程 (ARE)：

必须理解：

• Q, R 矩阵的物理含义（权重越大，对应的量越小）
• Riccati 方程的求解（`dlqr` / `care` 函数背后做了什么）
• LQR 是无限时域、无约束的；MPC 是有限时域、有约束的
• MPC 的终端代价 P 就是 LQR 的 Riccati 解

1.3 离散化方法

物理系统大多是连续的，但 MPC 在计算机上实现必须是离散的。

1.4 状态估计基础

实际 MPC 中状态往往不能直接测量，需要观测器。

• Luenberger 观测器（确定性系统）
• Kalman 滤波器（噪声系统）
• 理解观测器与 MPC 的分离原理

1.5 自查清单

进入第二层之前，确认你能：

• ✅ 写出任意线性系统的状态空间表示
• ✅ 手解离散 LQR 的 Riccati 方程（NumPy 实现）
• ✅ 用 ZOH 离散化一个连续系统
• ✅ 在 Simulink 或 Python 中仿真 LQR 闭环
• ✅ 描述 LQR 和 MPC 的核心区别

2 第二层：MPC 核心（3-4 周）

这是最关键的阶段。目标：手推 MPC 的 QP 形式 + 手写一个能跑的 MPC 控制器。

2.1 预测模型推导

给定系统 \(x_{k+1} = Ax_k + Bu_k\)，在时刻 \(k\) 已知 \(x_k = x_0\)：

写成矩阵形式：

其中：

这是 MPC 最重要的一组公式。必须亲手推一遍，理解每个矩阵块的来龙去脉。

2.2 QP 矩阵组装

把代价函数写成标准 QP 形式，这是 MPC 工程的精髓。

标准代价函数：

步骤1：定义块对角权重矩阵

步骤2：代入预测方程 \(\mathbf{x} = \mathcal{A}x_0 + \mathcal{B}\mathbf{u}\)

步骤3：得到标准 QP：

常数项 \(x_0^T \mathcal{A}^T \bar{Q} \mathcal{A} x_0\) 与优化变量无关，直接省略。

2.3 约束处理

这是 MPC 相对于 LQR 的核心优势。

输入约束：

状态约束（通过预测矩阵转化）：

变化率约束：

软约束：当硬约束导致 QP 无解时，引入松弛变量：

2.4 从零手写 MPC

这是最重要的环节。用 Python + OSQP 或 cvxpy 写一个能跑的 MPC。

import numpy as np
import cvxpy as cp

# 系统：双积分器 x = [p, v], u = a
dt = 0.1
A = np.array([[1, dt], [0, 1]])
B = np.array([[0.5*dt**2], [dt]])

# MPC 参数
N = 10  # 预测时域
Q = np.diag([10, 1])  # 状态权重
R = np.array([[0.1]]) # 控制权重
P = np.diag([100, 10]) # 终端权重

# 约束
u_max, u_min = 1, -1
x_max, x_min = np.array([5, 5]), np.array([-5, -5])

def solve_mpc(x0):
    x = cp.Variable((2, N+1))
    u = cp.Variable((1, N))

    cost = 0
    constraints = [x[:, 0] == x0]

    for i in range(N):
        cost += cp.quad_form(x[:, i], Q) + cp.quad_form(u[:, i], R)
        constraints += [x[:, i+1] == A @ x[:, i] + B @ u[:, i]]
        constraints += [u_min <= u[:, i], u[:, i] <= u_max]
        if i > 0:
            constraints += [x_min <= x[:, i], x[:, i] <= x_max]
    cost += cp.quad_form(x[:, N], P)
    constraints += [x_min <= x[:, N], x[:, N] <= x_max]

    prob = cp.Problem(cp.Minimize(cost), constraints)
    prob.solve(solver=cp.OSQP, warm_start=True)
    return u.value[:, 0], x.value, u.value

2.5 自查清单

• ✅ 手写出 \(\mathcal{A}\) 和 \(\mathcal{B}\) 预测矩阵
• ✅ 手算出 H 矩阵和 f 向量（从 Q, R, A, B, N 出发）
• ✅ 用 cvxpy / OSQP 实现无约束和有约束 MPC
• ✅ 能解释 N 太大/太小对系统的影响
• ✅ 能调 Q, R 达到期望的控制效果
• ✅ 理解软约束何时用、怎么用

3 第三层：NMPC 进阶（2-3 周）

当系统有强非线性时（倒立摆大角度摆动、机械臂、飞行器），线性 MPC 不够用，需要 NMPC。

3.1 非线性模型处理

NMPC 问题：

由于 \(f\) 是非线性的，QP 不再适用。需要求解非线性规划（NLP）。

三种主流方法：

3.2 NLP 求解算法

NMPC 的核心求解器：

3.3 从零手写 NMPC

import casadi as ca
import numpy as np

# 倒立摆非线性模型 (x = [theta, omega])
def pendulum_dynamics(x, u):
    g, L, m, b = 9.81, 1.0, 1.0, 0.1
    theta, omega = x[0], x[1]
    dtheta = omega
    domega = (g/L)*ca.sin(theta) + (1/(m*L**2))*u[0] - (b/m)*omega
    return ca.vertcat(dtheta, domega)

# RK4 离散化
def rk4_step(x, u, dt):
    k1 = pendulum_dynamics(x, u)
    k2 = pendulum_dynamics(x + dt/2*k1, u)
    k3 = pendulum_dynamics(x + dt/2*k2, u)
    k4 = pendulum_dynamics(x + dt*k3, u)
    return x + dt/6*(k1 + 2*k2 + 2*k3 + k4)

# 构建 NMPC
N, dt = 10, 0.05
x_sym = ca.MX.sym('x', 2)  # 符号变量
u_sym = ca.MX.sym('u', 1)
F = ca.Function('F', [x_sym, u_sym], [rk4_step(x_sym, u_sym, dt)])

opti = ca.Opti()  # NLP 问题
X = opti.variable(2, N+1)
U = opti.variable(1, N)

# 代价
Q, R = ca.diag([10, 1]), ca.diag([0.1])
cost = 0
for i in range(N):
    cost += X[:, i].T @ Q @ X[:, i] + U[:, i].T @ R @ U[:, i]
cost += X[:, N].T @ ca.diag([100, 10]) @ X[:, N]
opti.minimize(cost)

# 约束
opti.subject_to(X[:, 0] == [0.5, 0])  # 初始角度 0.5 rad
for i in range(N):
    opti.subject_to(X[:, i+1] == F(X[:, i], U[:, i]))
    opti.subject_to(U[:, i] <= 2)
    opti.subject_to(U[:, i] >= -2)

opti.solver('ipopt')
sol = opti.solve()
u_opt = sol.value(U[:, 0])

3.4 MPC vs NMPC 对比实验

同一个倒立摆系统，分别用线性 MPC（在平衡点线性化）和 NMPC 控制，对比效果：

3.5 自查清单

• ✅ 能用 CasADi 或 ACADOS 构建 NMPC
• ✅ 理解多重打靶法和直接配点法的区别
• ✅ 对比同一系统上的 MPC 和 NMPC 效果差异
• ✅ 了解 SQP 和内点法的基本迭代步骤
• ✅ 知道 RTI 为什么比标准 NMPC 快

4 第四层：深入理论（可选）

第四层是学术研究方向的进阶内容。如果你只是为了工程应用，第三层已经足够。

4.1 稳定性证明

MPC 的理论核心：证明闭环系统渐近稳定。

关键定理：如果满足以下条件，MPC 闭环渐近稳定：

1. 阶段代价 \(\ell(x,u) \geq \alpha(\|x\|)\)，正定（只在原点为零）
2. 终端代价 \(V_f\) 是终端局部控制器 \(\kappa_f(x)\) 的 Lyapunov 函数：
   $$V_f(f(x, \kappa_f(x))) - V_f(x) \leq -\ell(x, \kappa_f(x)), \quad \forall x \in \mathcal{X}_f$$
3. 终端约束集 \(\mathcal{X}_f\) 是 \(\kappa_f\) 的控制不变集
4. 初始时刻优化问题可行

证明思路：以最优代价 \(V_N^*(x_k)\) 作为 Lyapunov 候选函数，证明沿闭环轨迹单调下降：

这保证 \(x_k \to 0\)（渐近稳定）。

4.2 鲁棒 MPC

当系统存在模型不确定性或外部扰动时：

4.3 其他高级专题

📅 5. 学习计划与资源

5.1 建议时间线

5.2 推荐书目（按优先级）

5.3 代码资源

5.4 学习工具对比

❓ 6. 常见问题