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(k1) 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, i1) 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(i1, control_horizon)): Su[i*C.shape[0]:(i1)*C.shape[0], j*B.shape[1]:(j1)*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(PH, qf, AA, l-np.inf*np.ones(len(b)), ub, verboseFalse) 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, N1)) 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[:, k1] 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]) uRu return cost res minimize(nonlinear_mpc_cost, u_guess, args(x0, model, ref), boundsu_bounds, constraintsnonlinear_constraints)分布式MPC架构对于大规模系统可采用分解协调法将大系统分解为子系统分别优化博弈论方法各控制器作为博弈参与者达成均衡优先级分配按重要性分层优化鲁棒MPC设计考虑模型不确定性的方法包括最小-最大优化针对最坏情况设计随机MPC考虑概率分布Tube MPC设计不变集保证鲁棒性实际工程中MPC的性能高度依赖于模型精度。当模型存在显著误差时可结合数据驱动方法进行模型在线更新形成自适应MPC架构。