MPC模型预测控制实战：从理论到代码实现（Python示例）
MPC模型预测控制实战：从理论到代码实现（Python示例）
在工业控制和自动化领域，模型预测控制(MPC)已经成为处理多变量约束系统的主流方法。不同于传统的PID控制，MPC通过在线优化解决控制问题，特别适合处理具有时延、约束和交互影响的实际系统。本文将带您从零开始构建一个完整的MPC控制器，并通过Python代码展示如何将数学公式转化为可运行的算法。
1. MPC核心概念与实现框架
模型预测控制的魅力在于它将控制问题转化为一个滚动优化的数学问题。每次采样时刻，控制器都会：
- 基于当前状态和系统模型预测未来动态
- 求解一个有限时域的优化问题
- 仅应用优化解的第一个控制输入
- 在下一个采样时刻重复整个过程
这种"预测-优化-执行"的循环结构使MPC能够：
- 显式处理输入输出约束
- 自然地处理多变量耦合系统
- 通过重新优化补偿模型误差和干扰
典型MPC实现流程：
while not stop_condition: # 1. 获取当前状态 x = get_current_state() # 2. 求解优化问题 u_opt = solve_mpc_optimization(x) # 3. 应用第一个控制输入 apply_control(u_opt[0]) # 4. 等待下一个采样时刻 wait_for_next_sample()2. 系统建模与预测方程构建
任何MPC实现的第一步都是建立能够描述系统动态的数学模型。我们以离散状态空间模型为例：
$$ \begin{aligned} x(k+1) &= Ax(k) + Bu(k) \ y(k) &= Cx(k) \end{aligned} $$
其中$x\in\mathbb{R}^n$是状态向量，$u\in\mathbb{R}^m$是控制输入，$y\in\mathbb{R}^p$是系统输出。
预测方程推导步骤：
- 将模型转换为增量形式以减少稳态误差
- 递归应用状态方程构建多步预测
- 将预测方程整理为矩阵形式
def build_prediction_matrices(A, B, C, prediction_horizon, control_horizon): # 构建状态预测矩阵 Sx = np.vstack([np.linalg.matrix_power(A, i+1) for i in range(prediction_horizon)]) # 构建输入预测矩阵（下三角结构） Su = np.zeros((prediction_horizon*C.shape[0], control_horizon*B.shape[1])) for i in range(prediction_horizon): for j in range(min(i+1, control_horizon)): Su[i*C.shape[0]:(i+1)*C.shape[0], j*B.shape[1]:(j+1)*B.shape[1]] = C @ np.linalg.matrix_power(A, i-j) @ B return Sx, Su3. 优化问题构建与QP转化
MPC的核心是将控制问题转化为二次规划(QP)问题。典型的目标函数包含：
- 输出跟踪误差惩罚
- 控制量变化惩罚
目标函数数学表达：
$$ \min_{\Delta U} |Y - R|^2_Q + |\Delta U|^2_R $$
其中：
- $Y$是预测输出序列
- $R$是参考轨迹
- $\Delta U$是控制增量序列
- $Q$, $R$是权重矩阵
QP标准形式转化：
通过代数运算，可将目标函数转化为：
$$ \min_{\Delta U} \frac{1}{2}\Delta U^T H \Delta U + f^T \Delta U $$
其中：
- $H = 2(S_u^T Q S_u + R)$
- $f = -2S_u^T Q (R - S_x x_k)$
def mpc_to_qp(Sx, Su, Q, R, x_k, ref_traj): # 构建QP问题的H矩阵和f向量 H = 2 * (Su.T @ Q @ Su + R) f = -2 * Su.T @ Q @ (ref_traj - Sx @ x_k) return H, f4. 约束处理与求解器实现
实际系统中的物理限制必须作为约束条件加入优化问题。常见约束包括：
- 控制量幅值约束：$u_{min} \leq u \leq u_{max}$
- 控制增量约束：$\Delta u_{min} \leq \Delta u \leq \Delta u_{max}$
- 输出约束：$y_{min} \leq y \leq y_{max}$
约束矩阵构建示例：
def build_constraints(Su, u_prev, u_min, u_max, du_min, du_max, y_min, y_max): # 控制增量约束 A_du = np.vstack([np.eye(Su.shape[1]), -np.eye(Su.shape[1])]) b_du = np.hstack([np.tile(du_max, control_horizon), -np.tile(du_min, control_horizon)]) # 控制量幅值约束 L = np.tril(np.ones((control_horizon, control_horizon))) A_u = np.vstack([L, -L]) b_u = np.hstack([np.tile(u_max - u_prev, control_horizon), -np.tile(u_min - u_prev, control_horizon)]) # 输出约束 A_y = np.vstack([Su, -Su]) b_y = np.hstack([np.tile(y_max, prediction_horizon), -np.tile(y_min, prediction_horizon)]) # 合并所有约束 A = np.vstack([A_du, A_u, A_y]) b = np.hstack([b_du, b_u, b_y]) return A, bQP求解器选择：
Python中常用的QP求解器包括：
| 求解器 | 特点 | 适用场景 |
|---|---|---|
| CVXOPT | 精确但较慢 | 中小规模问题 |
| OSQP | 高效ADMM算法 | 大规模稀疏问题 |
| quadprog | 快速积极集法 | 中等规模问题 |
| SciPy | 内置但功能有限 | 简单问题原型 |
def solve_mpc(H, f, A, b): # 使用OSQP求解器示例 prob = osqp.OSQP() prob.setup(P=H, q=f, A=A, l=-np.inf*np.ones(len(b)), u=b, verbose=False) res = prob.solve() return res.x5. 完整MPC控制器实现与调参
将上述组件整合，我们得到完整的MPC控制器实现：
class MPCController: def __init__(self, A, B, C, pred_horizon, ctrl_horizon, Q, R): self.Sx, self.Su = build_prediction_matrices(A, B, C, pred_horizon, ctrl_horizon) self.Q = block_diag(*[Q]*pred_horizon) self.R = block_diag(*[R]*ctrl_horizon) self.pred_horizon = pred_horizon self.ctrl_horizon = ctrl_horizon self.u_prev = np.zeros(B.shape[1]) def compute_control(self, x, ref_traj, constraints): # 构建QP问题 H, f = mpc_to_qp(self.Sx, self.Su, self.Q, self.R, x, ref_traj) A, b = build_constraints(self.Su, self.u_prev, **constraints) # 求解QP du = solve_mpc(H, f, A, b) # 提取第一个控制增量并更新 du_opt = du[:self.Su.shape[1]//self.ctrl_horizon] u_opt = self.u_prev + du_opt self.u_prev = u_opt return u_opt关键调参建议：
预测时域与控制时域：
- 预测时域越长，控制性能越好但计算负担越重
- 控制时域通常为预测时域的1/3到1/2
权重矩阵选择：
- 输出误差权重(Q)决定跟踪精度
- 控制增量权重(R)影响控制平滑性
- 可通过Bryson规则初始化： $$ Q_{ii} = \frac{1}{y_{i,max}^2}, \quad R_{jj} = \frac{1}{u_{j,max}^2} $$
约束软化： 对关键输出约束添加松弛变量避免不可行问题：
# 在目标函数中添加松弛惩罚 H_aug = block_diag(H, 1e6*np.eye(num_constraints)) f_aug = np.hstack([f, np.zeros(num_constraints)])
6. 应用案例：水箱液位控制
考虑串联水箱系统的液位控制问题，其线性化模型为：
$$ \frac{d}{dt}\begin{bmatrix}h_1\h_2\end{bmatrix} = \begin{bmatrix}-\frac{1}{R_1A_1} & 0 \ \frac{1}{R_1A_2} & -\frac{1}{R_2A_2}\end{bmatrix} \begin{bmatrix}h_1\h_2\end{bmatrix} + \begin{bmatrix}\frac{1}{A_1}\0\end{bmatrix}q_{in} $$
离散化与MPC实现：
# 系统参数 A1, A2 = 0.5, 0.3 # 水箱截面积(m^2) R1, R2 = 0.2, 0.3 # 流阻(m^2/s) # 连续状态空间矩阵 Ac = np.array([[-1/(R1*A1), 0], [1/(R1*A2), -1/(R2*A2)]]) Bc = np.array([[1/A1], [0]]) Cc = np.eye(2) # 离散化(零阶保持) dt = 1.0 # 采样时间(s) A = expm(Ac*dt) B = np.linalg.inv(Ac) @ (A - np.eye(2)) @ Bc C = Cc # MPC参数 pred_horizon = 10 ctrl_horizon = 4 Q = np.diag([10, 5]) # 更关注h1的跟踪 R = np.array([[0.1]]) # 控制量变化惩罚 # 创建MPC控制器 mpc = MPCController(A, B, C, pred_horizon, ctrl_horizon, Q, R) # 模拟运行 N = 50 h = np.zeros((2, N+1)) u = np.zeros(N) ref = np.ones(2*pred_horizon) # 双水箱的参考轨迹 for k in range(N): # 获取当前状态 x = h[:, k] # MPC计算控制量 u[k] = mpc.compute_control(x, ref, { 'u_min': 0, 'u_max': 2, 'du_min': -0.5, 'du_max': 0.5, 'y_min': np.array([0, 0]), 'y_max': np.array([1.5, 1.5]) }) # 系统仿真 h[:, k+1] = A @ h[:, k] + B * u[k]性能优化技巧：
- 热启动：使用上一步的解作为当前优化的初始猜测
- 稀疏性利用：预测矩阵具有特定稀疏结构，可加速计算
- 实时性保障：
- 设置求解时间上限
- 准备备用控制策略（如PID）应对求解失败
7. 高级主题与扩展方向
非线性MPC处理：
- 连续线性化：在每个采样点重新线性化模型
- 直接非线性优化：使用IPOPT等求解器处理非线性问题
- 神经网络近似：用深度学习模型预测系统动态
# 非线性MPC示例框架 def nonlinear_mpc_cost(u_sequence, x0, model, ref): # 模拟系统响应 x = x0 cost = 0 for i, u in enumerate(u_sequence): x = nonlinear_model(x, u) # 非线性模型推进 cost += (x - ref[i])@Q@(x - ref[i]) + u@R@u return cost res = minimize(nonlinear_mpc_cost, u_guess, args=(x0, model, ref), bounds=u_bounds, constraints=nonlinear_constraints)分布式MPC架构：
对于大规模系统，可采用：
- 分解协调法：将大系统分解为子系统分别优化
- 博弈论方法：各控制器作为博弈参与者达成均衡
- 优先级分配：按重要性分层优化
鲁棒MPC设计：
考虑模型不确定性的方法包括：
- 最小-最大优化：针对最坏情况设计
- 随机MPC：考虑概率分布
- Tube MPC：设计不变集保证鲁棒性
实际工程中，MPC的性能高度依赖于模型精度。当模型存在显著误差时，可结合数据驱动方法进行模型在线更新，形成自适应MPC架构。