模型预测控制(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_steps1)] # 示例时间常数为5预测步长20 step_coeffs generate_step_coefficients(5, 20)对于多变量系统我们需要构建动态矩阵A其结构如下系数结构说明a₁第一时刻响应a₂第二时刻响应......aₚ第P时刻响应数学表达为y(k1) a₁Δu(k) y₀(k1) y(k2) a₁Δu(k1) a₂Δu(k) y₀(k2) ... y(kP) Σ[aᵢΔu(kP-i)] y₀(kP)3. 预测方程矩阵化实现将预测方程转化为矩阵形式可大幅提升计算效率特别适合Python的数值计算生态。我们构建的预测系统包含以下核心矩阵动态矩阵AP×M维包含阶跃响应系数控制增量ΔUM×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(i1, control_horizon)): A[i,j] step_coeffs[i-j] return A # 示例预测步长P5控制步长M3 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 const5. 带约束优化求解实际工程问题通常包含多种约束条件如控制量幅值限制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(k1) y(k1) - ŷ(k1|k) 校正 y_cor(ki) ŷ(ki) h_i·e(k1)实现代码示例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, P10, M5, q1, r0.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_max0.5, u_min-1, u_max1): 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--, labelReference) plt.plot(time_axis, outputs, b-, labelOutput) plt.ylabel(System Output) plt.legend() plt.subplot(2,1,2) plt.step(time_axis, controls, g-, wherepost, labelControl) plt.ylabel(Control Signal) plt.xlabel(Time (s)) plt.legend() plt.tight_layout() plt.show()关键调参经验预测步长P通常覆盖系统主要动态一般为上升时间的1.5-2倍控制步长M3-5步即可过大增加计算负担权重比q/r决定跟踪精度与控制成本的平衡约束设置应根据物理限制合理设置过紧会导致无解9. 工业应用中的实战技巧在实际工程部署MPC时有几个关键考量点模型失配处理定期更新阶跃响应系数增加扰动观测器放宽输出约束边界实时性保障# 预计算不变部分优化 H_inv np.linalg.inv(H) # 离线计算H的逆 def fast_solve(f): return -H_inv f.T # 在线仅需矩阵乘法数值稳定性增强对H矩阵进行正则化处理添加小量单位矩阵防止奇异采用带约束的QP求解器10. 进阶扩展方向对于更复杂的控制场景可考虑以下扩展非线性MPCfrom 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快变量高频控制慢变量低频更新分层优化架构分布式MPC子系统独立优化协调层处理耦合基于ADMM的分布式求解