简介面向需要理解和应用模型预测控制MPC的工程师与研究人员这份MATLAB实现资料包结合源码与文档完整覆盖MPC的核心流程建立系统动态模型、确定预测与控制时域、设定性能指标和约束条件、构建并求解优化问题以及滚动实施控制输入。压缩包共2个文件分别为一个可运行的MATLAB脚本和一份Word说明文档整体仅156KB内容集中且便于下载后快速阅读。说明文档对MPC原理、实现步骤及MATLAB相关函数进行了梳理包括系统模型构建、mpc控制器对象创建、模型仿真与结果分析等环节MATLAB脚本则提供算法实现框架读者可调整模型参数、预测时域或约束条件观察控制效果变化加深对滚动优化的理解。目前已有2196人学习浏览既适合MPC入门自学也可用于课程设计或课题预研是一份兼顾原理讲解与动手验证的实用资料。1. 模型预测控制是什么一个比PID多算一步的控制器模型预测控制算法实现听起来比PID要高级不少但核心逻辑其实就一句话每个采样周期都把未来一段时间的控制量算一遍只执行第一步到下一个周期再拿新测量值重算。我不止一次见到有人把MPC当成“调好Q和R就完事”的黑匣子结果一加约束就翻车。这篇文章用MATLAB作为落地载体从手写一个能跑的最小闭环开始到带输入约束和速度约束的双积分位置伺服实例再把预测时域、控制时域、采样周期这些参数讲透最后列出我实际调试里踩过的坑。适合想搞清楚MPC模型预测控制到底怎么落地而不是只停留在工具箱调用层面的工程师和学生。2. 核心三要素与MATLAB里能跑的最小MPC代码2.1 预测模型、滚动优化、反馈校正到底在算什么MPC模型预测控制在每个采样周期里只做三件事预测、优化、执行一步。预测模型负责回答“如果我这么控制未来Np拍系统会走到哪”滚动优化负责在约束条件下找一个让代价函数最小的控制序列反馈校正则是只把该序列的第一步作用到对象上下一拍再用新测量值重来。为什么一定要“滚动”因为预测模型永远不可能是对象的完美复制。模型失配、外部扰动、约束边界上的非线性都会让远期预测越来越不可信。滚动优化相当于每拍都把预测窗口重新往前推把误差及时修正回来。这也是MPC能处理约束的核心原因它在优化时就明确知道“未来输入不能超过1”而不是像PID那样等超调发生了再去拉回来。一个采样周期内的计算顺序是固定的读当前状态x(k) → 构建这一拍的二次规划问题 → 求解出最优控制序列 → 把序列第一个值u(k)送到执行器 → 下一拍重复。这里最关键的一点是QP的代价矩阵H通常可以提前算好每拍只需要更新线性项f这是后面性能优化的基础。2.2 用quadprog手写一个最小MPC目标函数与约束的矩阵构造先从一个一阶对象开始x(k1) a·x(k) b·u(k)跟踪常数设定值r输入约束|u| ≤ 1。这个对象足够简单能让你看清MPC每个矩阵在算什么而不是被状态空间的维度淹没。% 最小线性MPC: 一阶对象, 用quadprog做滚动优化 % 对象: x(k1) a*x(k) b*u(k) a 0.9; b 0.1; Np 5; Nc 3; % 预测时域5, 控制时域3 Q 1; R 0.5; % 构建预测矩阵 Phi(自由响应) 和 Gamma(输入响应) Phi zeros(Np,1); Gamma zeros(Np,Nc); for i 1:Np Phi(i) a^i; % 初值x(k)对未来第i拍的自由响应 for j 1:Nc if i j % 第j个控制增量对未来第i拍的响应, 注意是a^(i-j)*b Gamma(i,j) a^(i-j)*b; end end end % 代价函数 J (Y-r)*Q*(Y-r) u*R*u % quadprog求解的是 0.5*u*H*u f*u, 所以H和f都要按2倍折算 H 2*(Gamma*Q*Gamma R*eye(Nc)); Ain [eye(Nc); -eye(Nc)]; % |u|1 bin ones(2*Nc,1); x 0; r 2; u_rec zeros(50,1); x_rec zeros(50,1); for k 1:50 % 线性项里包含当前状态x和设定值r, 每拍都要更新 f -2*Gamma*Q*(r - Phi*x); u_opt quadprog(H, f, Ain, bin, [], [], -1, 1); u u_opt(1); % 只执行第一步 x a*x b*u; % 更新对象 u_rec(k) u; x_rec(k) x; end代码里最容易写错的是H和f的符号。quadprog默认求解min 0.5·u·H·u f·u而MPC代价展开后是u·(ΓQΓR)·u所以H必须乘2f必须带负号。如果你直接把ΓQΓ当作H传进去等于把权重缩小了一半闭环会明显变慢。另一个隐蔽的默认行为是Nc3意味着控制序列长度只有3但预测时域Np5那么第4、5拍的输入在预测模型里是0不是保持u(3)不变。很多从论文公式手写MPC的人在这里翻车——论文里常写“u(kj)u(kNc-1)”但MATLAB矩阵构造里不会自动这么干。需要保持时得把Gamma的后Np-Nc列补成“最后一步输入”的累积响应而不是空着。这个例子虽然小但已经具备完整MPC闭环预测矩阵、QP构造、滚动执行。跑完后检查x_rec应该稳定到2附近u_rec初始阶段会顶在1的边界上这就是约束在起作用。如果没有约束等效于一个带模型前馈的线性二次调节器这正好引出下一个问题为什么不用PID直接调。2.3 为什么先手写而不是直接上工具箱MATLAB的Model Predictive Control Toolbox里mpc对象配合mpcmove函数几行就能跑通一个MPC看起来省事得多。但我的建议是新手先手写一遍至少跑通上面这种最小例子。原因是工具箱把预测矩阵、QP变换、约束重整化都封装死了你只能看到Q、R、Np这几个旋钮一旦约束不可行或轨迹发散你根本不知道是哪个矩阵构造出了问题。手写一遍的价值在于你能亲眼看到自由响应Phi和输入响应Gamma是怎么叠加的速度约束是怎么变成Aineq·u ≤ bineq的H和f的符号约定为什么和quadprog差一个系数。这些一旦通了用工具箱只是换一个求解入口而已。实现方式适合场景代价手写quadprog学习原理、定制约束、部署到嵌入式要自己维护矩阵构造mpc对象标准线性MPC、快速原型验证黑匣子出问题难查nlmpc对象非线性模型、强非线性约束求解慢初值敏感工具箱适合在算法定型后做交叉验证这是我建议的工作顺序手写跑通 → 改约束和参数 → 工具箱对照 → 再回手写版本部署。下面一章就把这个流程套到一个真正有状态约束的实例上。3. 一个完整的模型预测控制实例双积分对象的位置伺服3.1 对象建模与离散化从连续模型到预测矩阵假设被控对象是小车在导轨上运动控制量是加速度u状态是位置x1和速度x2。连续状态方程写出来是dx/dt [0 1; 0 0]·x [0; 1]·u这是一个双积分对象名字很朴素但它是很多运动控制问题的标准抽象位置不能突变速度是位置的积分加速度直接由控制量决定。用MATLAB做离散化时我一般直接用c2d而不是手推欧拉公式因为零阶保持zoh离散化才是物理上D/A的真实行为。% 连续模型 - 零阶保持离散化 Ac [0 1; 0 0]; Bc [0; 1]; Cc [1 0]; Dc 0; Ts 0.1; sysd c2d(ss(Ac,Bc,Cc,Dc), Ts, zoh); Ad sysd.A; Bd sysd.B; % 结果是 Ad [1 Ts; 0 1], Bd [Ts^2/2; Ts]为什么不推荐欧拉离散在这个对象上欧拉会得到Bd [0; Ts]和真实零阶保持的[Ts^2/2; Ts]差了一个Ts/2的因子。Ts0.1时误差还不大但当你为了快速响应把Ts压到0.01欧拉的输入增益会偏小一半左右MPC的预测模型和实际对象对不上闭环就会带稳态误差甚至振荡。离散化之后预测方程就能写成一个统一的矩阵形式把未来Np步的状态全部表示成当前状态x(k)加控制序列U的函数。具体来说自由响应矩阵Phi的每一列是Ad^i·x(k)输入响应矩阵Gamma的每一列是Ad^(i-j)·Bd对第j个控制输入的响应。这个矩阵构造是后面所有约束和代价计算的地基。3.2 带输入和速度约束的MPC函数QP如何真正落地现在给这个实例加上两个硬约束加速度也就是输入u的绝对值不超过1速度x2的绝对值不超过0.5。这非常贴近实际——电机有最大出力小车有最大允许速度。MPC的亮点是它直接在优化问题里声明这些硬约束而不是像PID那样把限幅放在控制器外面事后截断。function [U_opt, feas] mpc_pos_servo(xk, r, params) % 模型预测控制实例: 双积分位置伺服 % 状态 x[位置;速度], 输入 u 为加速度 Ad params.Ad; Bd params.Bd; Np params.Np; Nc params.Nc; Q params.Q; R params.R; % 预测矩阵: Phi(2Np x 2), Gamma(2Np x Nc) Phi zeros(2*Np, 2); Gamma zeros(2*Np, Nc); for i 1:Np Phi(2*i-1:2*i, :) Ad^i; for j 1:Nc if j i Gamma(2*i-1:2*i, j) Ad^(i-j)*Bd; end end end % 只看位置分量: 用矩阵E把奇数行抽出来 E zeros(Np, 2*Np); for i 1:Np E(i, 2*i-1) 1; end Gy E*Gamma; Phiy E*Phi*xk; % 代价和quadprog标准形式 H 2*(Gy*Q*Gy R*eye(Nc)); f -2*Gy*Q*(r - Phiy); % 输入约束 |u|1, 作用在Nc个控制量上 A_u [eye(Nc); -eye(Nc)]; b_u ones(2*Nc,1); % 速度约束 |v|0.5, 作用在预测域每一拍上 % 先构造Sv, 把每个预测步的速度分量取出来 Sv zeros(Np, 2*Np); for i 1:Np Sv(i, 2*i) 1; end A_v [Sv*Gamma; -Sv*Gamma]; b_v [0.5 - Sv*Phi*xk; 0.5 Sv*Phi*xk]; % 合并所有不等式约束, 上下界给输入限幅 A [A_u; A_v]; b [b_u; b_v]; opts optimoptions(quadprog, Display, off); [U_opt, ~, exitflag] quadprog(... H, f, A, b, [], [], -ones(Nc,1), ones(Nc,1), [], opts); feas (exitflag 1); end这里最需要理解的是速度约束那一坨矩阵。Sv·Gamma·U是控制序列对每个预测步速度的贡献Sv·Phi·xk是当前状态对每个预测步速度的自由贡献。约束写成 Sv·Gamma·U ≤ 0.5 - Sv·Phi·xk 和 -Sv·Gamma·U ≤ 0.5 Sv·Phi·xk合到一起就是A_v·U ≤ b_v。注意b_v里包含状态当前值所以这一步每拍都得重新算不能提前缓存。输入约束里的上下界[-1, 1]其实是可选的因为Ain已经包含了|u|≤1。但加上它能让quadprog更快收敛也能防止数值误差导致边界穿透。exitflag用来判断这一拍到底有没有找到可行解这直接关系到后面主仿真脚本的兜底策略。3.3 闭环仿真脚本先算最优再更新状态有了MPC求解函数主循环就非常清爽了。每次先根据当前状态出最优控制序列取第一个值作用到对象上然后用离散模型更新状态。注意顺序一定是先求u再更新x反过来就等于把最优控制用在了旧状态上这在快速系统上会带来一个采样周期的滞后。% 主脚本: 闭环仿真 params.Ad Ad; params.Bd Bd; params.Np 20; params.Nc 4; params.Q 10; params.R 1; x [0; 0]; r 2; % 从原点出发, 跟踪到位置2 Nsim 80; x_log zeros(Nsim,2); u_log zeros(Nsim,1); for k 1:Nsim [u_opt, ok] mpc_pos_servo(x, r, params); if ~ok % 这一拍不可行, 用上一拍控制量兜底 warning(第%d拍QP不可行, 回退到上一拍输入, k); u_opt u_log(max(k-1,1)) * ones(params.Nc,1); end u_apply u_opt(1); % 只执行第一步 x Ad*x Bd*u_apply; % 更新对象状态 x_log(k,:) x; u_log(k) u_apply; end t (1:Nsim)*Ts; plot(t, x_log(:,1), b-, LineWidth, 1.5); hold on plot(t, u_log, r-, LineWidth, 1); legend(位置, 控制量); xlabel(时间(s)); grid on跑完看三件事。第一位置是否收敛到2有没有明显的超调或等幅振荡。第二u_log初始阶段是不是顶在1的边界上——如果顶住了再松开说明输入约束确实被优化器主动处理了而不是截断。第三把速度曲线也画出来看它在0.5附近是否被“削平”这是MPC区别于PID限幅的最直观证据速度约束是在预测域内提前预判并减速的不是到了界限才生硬地夹断。如果你把Nc改成1再跑会发现响应明显变“肉”因为优化器只决定当前一步后面全靠自由响应。把Np从20改成5约束激活的位置会往后推动态变得迟钝。这些现象都是下一章参数调试的引子。4. MPC参数怎么调预测时域、控制时域与权重的落地经验4.1 预测时域Np太长太短分别会怎样预测时域Np决定MPC向前看多远单位是采样周期。Np太短最典型的症状是约束被“无视”——优化器只看到眼前几步速度约束还没来得及起效下一拍状态已经冲出去了。我见过有人把Np设成3闭环仿真里速度反复撞上约束边界看起来像极限环实际上不是就是预测窗口太短。Np太长的副作用有两个。一是QP变量不变但约束行数和预测矩阵变大计算量上升二是远期预测在模型失配时根本不靠谱你把权重压在Np最后几步上等于在为一个注定不准的遥远未来做决策动态响应反而偏慢。经验取值让Np·Ts覆盖系统开环主导时间常数的3到5倍。对双积分这种没有天然时间常数的对象就从“从0加速到设定值需要多少个采样周期”估算然后乘2到3。本文实例里Ts0.1期望上升时间约1秒Np取20正好覆盖2秒这就是比较稳的起点。4.2 控制时域Nc与权重Q、R的调整顺序控制时域Nc是优化变量里“自由”控制量的个数。Nc1时控制器最保守因为它只决定当前一步后面全部按0处理。Nc越大控制器越激进但Hessian矩阵变大求解时间上升而且对模型失配更敏感。我在运动控制里常用的起点是Np的15%到25%也就是Np20时取Nc3到5。权重Q和R的调整顺序比绝对值更重要。我的习惯是先固定R1把Q从小到大扫一遍Q很小时系统懒洋洋Q增大到临界点后开始出现输入饱和或振荡然后再回头放大R、减小Q/R比值。切忌一上来Q给100、R给0.01那样MPC会像高增益PID一样疯狂动作约束根本压不住。参数作用经验起点调大后的表现Np预测窗口长度时间常数的3~5倍计算变慢、响应变钝Nc自由控制量个数Np的15%~25%更激进但易抖Q输出跟踪权重先固定为1跟踪更快、输入更猛R输入惩罚权重先固定为1输入更平滑、响应更慢Ts采样周期带宽的10~20倍数值条件差或约束穿透一个很多人不知道的小技巧权重只对“趋势”敏感不要指望靠微调Q/R把振荡去掉。振荡往往是Nc太大或Np太短权重反而背了锅。4.3 采样周期对可行域的影响采样周期Ts在MPC里不只是离散化步长它直接决定预测矩阵的条件数和约束的“可见度”。Ts太小Bd项变成Ts²/2的量级Hessian里R这一项相对过大quadprog容易报数值病态同时覆盖同样物理时段需要的Np大幅增加求解时间变长。Ts太大约束只能在采样点上被检查两个采样点之间速度可能已经超限这就违背了“约束满足”的初衷。我的经验是Ts取被控系统期望带宽的10到20倍采样率。对一个期望闭环带宽2 rad/s的系统Ts取0.02到0.1都合理再小就是为计算量买单。在本文双积分实例里Ts0.1配合Np20已经能让速度约束在预测域内被提前看到如果Ts变成0.5Np就得拉到10以上才能覆盖同样时间否则约束刚看到就撞上了。采样周期还影响一个容易被忽略的点执行器的零阶保持行为。MPC算出的u是采样周期内的恒定值如果执行器带宽远低于1/Ts实际输出跟不上一拍一个值的指令相当于在对象前面串了一个额外滞后。这时不要把Ts继续调小反而应该放大Ts或把执行器动态写进预测模型。5. 模型预测控制实现中的避坑与常见问题排查5.1 现象闭环仿真前几步就发散发NaN原因最常见的是离散化方式不一致。有人用c2d得到Ad、Bd但仿真更新时又手写欧拉公式或者预测模型用欧拉、对象模型用零阶保持两边增益对不上MPC算出的控制量对真实对象过度激励。另一个常见原因是quadprog返回NaN却没用exitflag检查。解决统一用c2d(zoh)得到的Ad、Bd同时做预测和仿真。主循环里加一行feas检查exitflag不等于1就打印状态和约束上界不要闷头跑。5.2 现象quadprog报“No feasible solution found”原因硬约束本身就不相容。最典型的是初始速度已经在0.5以上下一拍还要等加速度把速度拉回来但预测域内每一拍都强制|v|≤0.5第一拍就不可行。另一个是把输入约束和状态约束同时收紧优化变量根本没有空地。解决先去掉速度约束跑通确认MPC本身没问题再把速度约束只加在预测域的后半段比如i≥3才开始生效最后给速度约束加软约束也就是在代价函数里引入一个slack变量并惩罚它让约束“尽量满足”而不是“必须满足”。软约束的代价是用户要调惩罚系数rhorho太小约束会明显越界太大又回到硬约束。5.3 现象模型失配下稳态误差一直消不掉原因MPC本质是内模控制预测模型增益偏了闭环稳态输出就会偏。比如真实对象增益比模型大10%MPC算出的输入对应的输出就偏低10%你没有PID里的积分项去消除这个余差。解决两个办法。一是在预测模型里加不可测扰动估计用当前测量值y(k)和上一拍预测值y(k|k-1)的差作为扰动修正未来Np步的预测输出这是MPC里的标准反馈校正。二是改用增量型MPC把优化变量从u(k)换成Δu(k)代价函数里惩罚Δu等于引入了积分作用。增量型改造要增广状态把u(k-1)加进状态向量矩阵构造会复杂一些但工程上非常值得。5.4 现象MATLAB里跑得很慢一帧要几十毫秒原因很多人把H、A、b这些不随状态变化的矩阵放在MPC函数里每拍重建一遍quadprog每次都要重新分配内存。还有一个隐藏问题主循环里用了动态增长的数组MATLAB在循环里不停复制数组。解决把Phi、Gamma、E、Sv、A_u、b_u全部提到函数外面作为结构体缓存只把H和约束矩阵中与x(k)相关的部分在函数里更新。仿真循环前用zeros预分配x_log、u_log。这些改完同样算例耗时通常能降到原来的十分之一。如果还要更快就把优化问题转化成mex或MATLAB Coder部署但那属于另一个话题。5.5 现象nlmpc仿真卡死或迭代极慢原因非线性MPC的求解器对初值极敏感。默认初值离可行域太远或者优化器在每一步都从零开始迭代几十次还没靠近最优解。另一个常见原因是状态估计器初值与真实状态差太大预测模型一路错下去。解决给nlmpc指定InitialGuess或设置InitialState用上一拍的最优解作为这一拍的初始猜测先用无约束线性MPC的结果跑一遍把轨迹存下来作为非线性MPC的热启动初值最后把OptimizationOptions里的MaxIterations放宽到200但StepTolerance别设太小否则会陷入振荡迭代。6. 进阶终端代价与MATLAB工具箱的交叉验证6.1 用mpc工具箱快速对比手写结果当你手写的MPC能稳定跑通双积分实例后强烈建议做一次交叉验证用Model Predictive Control Toolbox的mpc对象搭同样的预测时域、控制时域和权重对比两条闭环轨迹。如果形状不一致说明手写代码的某个矩阵构造有偏差如果完全重合你已经具备脱离工具箱独立实现MPC的能力。% 用MPC工具箱搭建同样的双积分对象 mpcobj mpc(ss(Ad,Bd,Cc,Dc), Ts); mpcobj.PredictionHorizon 20; mpcobj.ControlHorizon 4; mpcobj.Weights.OutputVariables 10; mpcobj.Weights.ManipulatedVariables 1; % 给控制量和状态加约束 mpcobj.ManipulatedVariables.Min -1; mpcobj.ManipulatedVariables.Max 1; sim(mpcobj, Nsim, r);工具箱的作用不只是验证。当你后面要做模型失配分析、状态估计器设计或C代码生成时mpcobj提供了一套标准化接口比手写矩阵鲁棒得多。但它始终是个黑匣子一旦遇到“约束不可行”或“软约束穿透”这类问题还是要回手写版本里看每个约束矩阵的数值所以两条线都要留着。6.2 一个具体技巧把终端代价换成离散Lyapunov解最后一个值得做的改进是终端代价。标准的有限时域MPC只保证“在窗口内最优”不保证闭环稳定性。给预测末端的终端状态加一个权重矩阵P相当于告诉优化器“结束那一刻的状态也需要被惩罚”这能让稳定性从工程调参变成有理论依据。P的推荐取法是解离散Riccati方程MATLAB的dlqr可以直接出。% 终端代价: 用dlqr解出终端权P, 把它折算进H和f Qc 10; Rc 1; % 终端代价的设计权重 [~, P, ~] dlqr(Ad, Bd, Qc, Rc); % 终端状态 x(kNp) Phi_end*xk Gamma_end*U Phi_end Ad^Np; Gamma_end zeros(1, Nc); for j 1:Nc Gamma_end(j) Ad^(Np-j)*Bd(1); % 位置端部响应 end % 把 Phi_end*P*Phi_end 和 Gamma_end*P*Gamma_end 加入H和f H H 2*(Gamma_end*P*Gamma_end); f f - 2*(r - Phi_end*xk)*P*Gamma_end;这段代码的核心思想是把“最后一步状态离目标有多远”直接放进代价让优化器在窗口尽头依然有收敛意图。调P的时候注意Qc、Rc是终端代价自己的权重不需要和主代价的Q、R保持一致它们影响的是“收尾的力度”。我做第一个MPC时习惯把Q调到很大觉得跟踪快就是好结果输入饱和、约束穿透一起出现折腾了两天才意识到可行域才是瓶颈。后来我养成了固定R、从小到大扫Q的习惯每改一组参数先看约束激活点再下结论。这个方向最大的门槛不是公式推导而是你愿不愿意把每一个矩阵、每一个符号的约定都亲手验证一遍。希望帮到你。本文还有配套的精品资源点击获取