用MATLAB手把手实现一个简易的MPC控制器附完整代码与避坑指南在工业控制和自动化领域模型预测控制MPC因其出色的多变量处理能力和约束处理特性而备受青睐。但对于初学者来说从理论到实践的跨越往往充满挑战——复杂的矩阵运算、晦涩的优化问题还有那些让人头疼的维度匹配错误。本文将用最直观的方式带你用MATLAB从零搭建一个可运行的MPC控制器避开那些教科书不会告诉你的坑。1. 基础准备理解MPC的核心骨架MPC本质上是一个预测-优化-执行的循环过程。想象你在驾驶汽车时会不断预测前方路况调整方向盘和油门这就是MPC的直观体现。数学上这个循环包含三个关键环节预测模型用状态空间方程描述系统动态滚动优化在每个时刻求解最优控制序列反馈校正用实际测量值修正预测误差对于离散线性系统状态空间方程可以表示为% 离散状态空间方程示例 Ad [0.8 0.2; -0.1 0.9]; % 状态矩阵 Bd [0.5; 0.3]; % 输入矩阵 Cd [1 0]; % 输出矩阵提示实际工程中通常先用c2d函数将连续模型离散化采样时间的选择会影响控制性能2. 预测矩阵构建从数学到代码预测时域内的系统行为可以表示为X F * x0 Φ * U E其中F和Φ是最关键的预测矩阵。让我们用MATLAB构建它们function [F, Phi, E] build_prediction_matrices(Ad, Bd, Cd, Dd, Np, Nc) nx size(Ad,1); % 状态维度 nu size(Bd,2); % 输入维度 ny size(Cd,1); % 输出维度 F zeros(Np*ny, nx); Phi zeros(Np*ny, Nc*nu); E zeros(Np*ny, size(Dd,2)); % 构建F矩阵 F(1:ny,:) Cd * Ad; for k 2:Np F((k-1)*ny1:k*ny,:) F((k-2)*ny1:(k-1)*ny,:) * Ad; end % 构建Φ矩阵 for i 1:Np for j 1:min(i,Nc) Phi((i-1)*ny1:i*ny, (j-1)*nu1:j*nu) ... Cd * Ad^(i-j) * Bd; end end % 构建扰动矩阵E假设有常数扰动 E(1:ny,:) Cd * Dd; for k 2:Np E((k-1)*ny1:k*ny,:) Cd * Ad^(k-1) * Dd E((k-2)*ny1:(k-1)*ny,:); end end常见错误矩阵维度不匹配通常源于Ad,Bd,Cd维度定义不一致建议在代码开头添加维度检查3. 优化问题求解QP问题的实战转换MPC的核心是将控制问题转化为二次规划(QP)min J U*H*U 2*U*G s.t. U_min ≤ U ≤ U_max对应的MATLAB实现function [U_opt, flag] solve_mpc_qp(H, G, U_min, U_max) nu length(U_min); Nc nu / size(U_min,1); % 构建不等式约束 A*U ≤ b A [eye(nu); -eye(nu)]; b [U_max; -U_min]; % 使用quadprog求解 options optimoptions(quadprog, Display, none); [U_opt, ~, exitflag] quadprog(H, G, A, b, [], [], [], [], [], options); flag exitflag 1; if ~flag warning(QP求解失败尝试无约束解); U_opt -H\G; % 退化为无约束解 end end参数选择经验值参数作用典型取值Np预测时域10-20Nc控制时域3-5Q状态权重对角矩阵(1-10)R输入权重小值(0.1-1)4. 闭环实现与调试技巧完整的MPC闭环实现需要处理几个关键点function [u, x_pred] mpc_controller(x, u_prev, Ad, Bd, Cd, Q, R, Np, Nc, umax, umin) % 构建预测矩阵 [F, Phi, E] build_prediction_matrices(Ad, Bd, Cd, 0, Np, Nc); % 构建QP参数 Qe kron(eye(Np), Q); Re kron(eye(Nc), R); H Phi*Qe*Phi Re; G Phi*Qe*(F*x); % 求解QP U_min kron(ones(Nc,1), umin); U_max kron(ones(Nc,1), umax); [U_opt, ~] solve_mpc_qp(H, G, U_min, U_max); % 取第一个控制量 u U_opt(1:length(umin)); x_pred F*x Phi*U_opt; end调试中常见问题及解决方案QP求解失败检查H矩阵的正定性尝试减小Nc或放松约束控制发散增加Q权重强化状态调节检查模型(Ad,Bd)是否正确计算延迟减少Np/Nc使用热启动技巧用上一时刻的解作为初始值% 典型的主循环结构 for k 1:sim_steps [u, x_pred] mpc_controller(x, u, Ad, Bd, Cd, Q, R, Np, Nc, umax, umin); x Ad*x Bd*u; % 系统动态 y Cd*x; % 记录数据 U_history(:,k) u; X_history(:,k) x; end5. 进阶技巧提升MPC性能的实用方法当基础MPC运行稳定后可以考虑以下优化1. 终端代价与约束% 在H矩阵中添加终端代价 P dare(Ad, Bd, Q, R); % 求解Riccati方程 H Phi*Qe*Phi Re Phi(end-ny1:end,:)*P*Phi(end-ny1:end,:);2. 扰动估计与抗偏移% 增广状态空间包含扰动估计 Ad_aug [Ad Bd; zeros(nu,nx) eye(nu)]; Bd_aug [Bd; eye(nu)];3. 软约束处理% 在QP中引入松弛变量 H blkdiag(Phi*Qe*Phi Re, 1e6*eye(ny)); G [Phi*Qe*F*x; zeros(ny,1)]; A [A zeros(size(A,1),ny); -Cd*Phi -eye(ny)]; b [b; -ymin Cd*F*x];6. 完整代码实现与测试案例以下是一个完整的SISO系统MPC控制示例%% 系统定义 sys ss(tf([1],[1 1.5 1])); % 二阶振荡系统 Ts 0.1; % 采样时间 sysd c2d(sys, Ts); % 离散化 Ad sysd.A; Bd sysd.B; Cd sysd.C; %% MPC参数 Np 15; % 预测时域 Nc 3; % 控制时域 Q 10; % 输出权重 R 0.1; % 输入权重 umax 1; umin -1; % 输入约束 %% 仿真设置 sim_time 5; % 秒 sim_steps sim_time/Ts; ref 1; % 参考值 %% 初始化 x zeros(size(Ad,1),1); u 0; X_history zeros(size(Ad,1), sim_steps); U_history zeros(1, sim_steps); Y_history zeros(1, sim_steps); %% 主循环 for k 1:sim_steps % 调用MPC控制器 [u, ~] mpc_controller(x, u, Ad, Bd, Cd, Q, R, Np, Nc, umax, umin); % 系统动态 x Ad*x Bd*u; y Cd*x; % 记录数据 U_history(:,k) u; X_history(:,k) x; Y_history(:,k) y; end %% 可视化 t (0:sim_steps-1)*Ts; figure; subplot(2,1,1); plot(t,Y_history); ylabel(Output); subplot(2,1,2); plot(t,U_history); ylabel(Control); xlabel(Time(s));这个案例展示了如何控制一个二阶振荡系统跟踪阶跃信号。在实际项目中你可能需要添加参考跟踪项实现状态估计器如Kalman滤波处理输入速率约束加入积分动作消除稳态误差7. 性能优化与实时实现当需要在嵌入式设备上部署MPC时考虑以下优化策略代码生成友好实现function [u, status] mpc_autocode(x, u_prev, params) % 将矩阵计算展开为标量运算 % 适用于代码生成 persistent H cholH % 复用矩阵分解 if isempty(H) % 离线计算不变部分 [~, Phi, ~] build_prediction_matrices(... params.Ad, params.Bd, params.Cd, 0, params.Np, params.Nc); H Phi*kron(eye(params.Np),params.Q)*Phi kron(eye(params.Nc),params.R); cholH chol(H,lower); end % 在线计算 Fx compute_Fx(x, params); % 展开的矩阵乘法 G Phi*kron(eye(params.Np),params.Q)*Fx; % 求解三角系统 z linsolve(cholH, G, struct(LT,true)); U_opt -linsolve(cholH, z, struct(UT,true)); u U_opt(1) params.Kf*u_prev; % 可能添加滤波 status 1; end计算耗时优化技巧预计算不变矩阵如H矩阵的Cholesky分解使用定点算术减少计算量采用显式MPC预先计算最优控制律的分区表达式在实时性要求高的场合实测表明对于4状态2输入的系统在STM32H7上单次QP求解可控制在1ms以内。