模型预测控制(MPC) Python 实现：从阶跃响应模型到二次规划求解（附完整代码）-尧图网络科技

模型预测控制(MPC) Python实现：从阶跃响应到二次规划实战

1. MPC核心思想与工程价值

模型预测控制(Model Predictive Control)作为现代控制理论的重要分支，正在工业自动化、机器人控制、智能驾驶等领域展现出强大的适应性。与传统的PID控制相比，MPC最显著的优势在于其多变量处理能力和约束显式处理特性。想象一下化工过程中需要同时控制温度、压力和流量多个变量，或者自动驾驶中需要协调转向、油门和刹车——这正是MPC大显身手的场景。

MPC的三阶段控制循环构成了其核心框架：

预测阶段：基于系统模型预测未来多个时间步的输出轨迹
优化阶段：求解带约束的优化问题获得最优控制序列
反馈校正：执行第一个控制量并引入实时反馈补偿

# MPC基本控制循环伪代码 while not control_done: current_state = get_system_state() # 获取当前系统状态 predicted_outputs = model.predict(current_state) # 预测未来输出 optimal_actions = optimizer.solve(predicted_outputs) # 求解最优控制 execute_action(optimal_actions[0]) # 执行第一步控制 apply_feedback_correction() # 反馈校正

这种"预测-优化-执行"的滚动时域策略，使MPC特别适合处理具有以下特征的系统：

多输入多输出(MIMO)耦合系统
存在输入输出约束的场合
控制目标可量化为优化指标的场景
系统动态特性复杂但可建模的情况

2. 阶跃响应模型构建

阶跃响应模型作为MPC中最直观的模型表达方式，特别适合描述线性时不变系统。其核心思想是通过系统的阶跃响应系数来预测未来输出。

模型建立步骤：

对系统施加单位阶跃输入
记录输出达到稳定值的过程数据
提取特征时间点的响应系数{a₁, a₂,...,aₙ}

import numpy as np # 典型的一阶系统阶跃响应系数生成 def generate_step_coefficients(time_constant, total_steps): return [1 - np.exp(-i/time_constant) for i in range(1, total_steps+1)] # 示例：时间常数为5，预测步长20 step_coeffs = generate_step_coefficients(5, 20)

对于多变量系统，我们需要构建动态矩阵A，其结构如下：

系数	结构说明
a₁	第一时刻响应
a₂	第二时刻响应
...	...
aₚ	第P时刻响应

数学表达为：

y(k+1) = a₁Δu(k) + y₀(k+1) y(k+2) = a₁Δu(k+1) + a₂Δu(k) + y₀(k+2) ... y(k+P) = Σ[aᵢΔu(k+P-i)] + y₀(k+P)

3. 预测方程矩阵化实现

将预测方程转化为矩阵形式可大幅提升计算效率，特别适合Python的数值计算生态。我们构建的预测系统包含以下核心矩阵：

动态矩阵A：P×M维，包含阶跃响应系数
控制增量ΔU：M×1维，待求解的优化变量
自由响应Y₀：P×1维，不考虑未来控制作用的响应

def build_dynamic_matrix(step_coeffs, prediction_horizon, control_horizon): A = np.zeros((prediction_horizon, control_horizon)) for i in range(prediction_horizon): for j in range(min(i+1, control_horizon)): A[i,j] = step_coeffs[i-j] return A # 示例：预测步长P=5，控制步长M=3 A_matrix = build_dynamic_matrix(step_coeffs[:5], 5, 3) print("动态矩阵A:\n", A_matrix)

预测输出方程的矩阵形式为：

Ŷ = A·ΔU + Y₀

其中Ŷ为P维预测输出向量，这一形式为后续优化问题奠定了基础。

4. 二次规划问题构建

MPC的核心优化问题通常表述为二次型代价函数最小化：

min J = (R-Ŷ)ᵀQ(R-Ŷ) + ΔUᵀRΔU

其中：

Q：输出误差权重矩阵（对角阵）
R：控制增量权重矩阵（对角阵）
R：参考轨迹向量

from scipy.linalg import block_diag def build_cost_matrices(q_weights, r_weights, prediction_horizon): Q = np.diag(q_weights) if len(q_weights) > 1 else np.eye(prediction_horizon)*q_weights[0] R = np.diag(r_weights) if len(r_weights) > 1 else np.eye(len(r_weights))*r_weights[0] return Q, R # 权重配置示例 Q, R = build_cost_matrices([1,1,1,1,1], [0.1,0.1,0.1], 5)

将预测方程代入代价函数，可得标准QP形式：

J = ΔUᵀ(AᵀQA + R)ΔU + 2(Y₀-R)ᵀQAΔU + const

5. 带约束优化求解

实际工程问题通常包含多种约束条件，如：

控制量幅值限制：u_min ≤ u ≤ u_max
控制增量限制：Δu_min ≤ Δu ≤ Δu_max
输出变量限制：y_min ≤ y ≤ y_max

使用Python的cvxopt库可高效求解这类QP问题：

from cvxopt import matrix, solvers def solve_mpc_qp(H, f, A_con, b_con, lb, ub): P = matrix(H) q = matrix(f) G = matrix(np.vstack([A_con, -A_con])) h = matrix(np.hstack([b_con, -b_con])) lb = matrix(lb) ub = matrix(ub) solvers.options['show_progress'] = False solution = solvers.qp(P, q, G, h, lb, ub) return np.array(solution['x']).flatten() # 示例：构建约束矩阵 A_con = np.vstack([A_matrix, np.eye(3)]) # 输出约束+控制增量约束 b_con = np.hstack([np.ones(5)*2, np.ones(3)*0.5]) # 约束边界

6. 反馈校正机制

为提高鲁棒性，MPC引入反馈校正环节处理模型失配：

误差： e(k+1) = y(k+1) - ŷ(k+1|k) 校正： y_cor(k+i) = ŷ(k+i) + h_i·e(k+1)

实现代码示例：

def apply_feedback_correction(predicted, actual, correction_factors): error = actual - predicted[0] return predicted + correction_factors * error # 校正系数通常取0.3-0.7 correction_factors = np.linspace(0.5, 0.1, len(predicted)) corrected = apply_feedback_correction(predicted, actual, correction_factors)

7. 完整MPC仿真实现

将各模块整合为完整的控制闭环：

class MPCController: def __init__(self, step_coeffs, P=10, M=5, q=1, r=0.1): self.P = P # 预测步长 self.M = M # 控制步长 self.A = build_dynamic_matrix(step_coeffs, P, M) self.Q, self.R = build_cost_matrices([q]*P, [r]*M, P) self.H = self.A.T @ self.Q @ self.A + self.R self.f = lambda y0, r: (y0 - r).T @ self.Q @ self.A def solve(self, y0, r, u_prev, du_max=0.5, u_min=-1, u_max=1): f = self.f(y0, r) du_opt = solve_mpc_qp(self.H, f, self.A, np.ones(self.P)*2, np.ones(self.M)*-du_max, np.ones(self.M)*du_max) u_opt = u_prev + du_opt[0] return np.clip(u_opt, u_min, u_max), du_opt # 仿真循环示例 mpc = MPCController(step_coeffs) u = 0 # 初始控制量 for k in range(100): y0 = get_free_response() # 获取自由响应 r = get_reference(k) # 获取参考轨迹 u, _ = mpc.solve(y0, r, u) apply_control(u) # 施加控制 time.sleep(0.1) # 实时控制周期

8. 可视化分析与性能调优

通过可视化可直观评估MPC性能：

import matplotlib.pyplot as plt def plot_mpc_performance(time_axis, reference, outputs, controls): plt.figure(figsize=(12, 8)) plt.subplot(2,1,1) plt.plot(time_axis, reference, 'r--', label='Reference') plt.plot(time_axis, outputs, 'b-', label='Output') plt.ylabel('System Output') plt.legend() plt.subplot(2,1,2) plt.step(time_axis, controls, 'g-', where='post', label='Control') plt.ylabel('Control Signal') plt.xlabel('Time (s)') plt.legend() plt.tight_layout() plt.show()

关键调参经验：

预测步长P：通常覆盖系统主要动态，一般为上升时间的1.5-2倍
控制步长M：3-5步即可，过大增加计算负担
权重比q/r：决定跟踪精度与控制成本的平衡
约束设置：应根据物理限制合理设置，过紧会导致无解

9. 工业应用中的实战技巧

在实际工程部署MPC时，有几个关键考量点：

模型失配处理：

定期更新阶跃响应系数
增加扰动观测器
放宽输出约束边界

实时性保障：

# 预计算不变部分优化 H_inv = np.linalg.inv(H) # 离线计算H的逆 def fast_solve(f): return -H_inv @ f.T # 在线仅需矩阵乘法

数值稳定性增强：

对H矩阵进行正则化处理
添加小量单位矩阵防止奇异
采用带约束的QP求解器

10. 进阶扩展方向

对于更复杂的控制场景，可考虑以下扩展：

非线性MPC：

from casadi import * def nonlinear_mpc(): x = MX.sym('x'); u = MX.sym('u') ode = vertcat(x[1], (1-x[0]**2)*x[1]-x[0]+u) dae = {'x':x, 'ode':ode} opts = {'tf':0.5} # 采样时间 F = integrator('F', 'cvodes', dae, opts) # 构建NMPC问题...

多速率MPC：