简介《模型预测控制理论、计算与设计》第二版是MPC领域公认的经典教材面向控制工程专业学生、科研人员与工程师系统串联预测控制从理论推导到算法实现再到控制器设计的完整链路。第二版新增第8章「数值最优控制」由Moritz M. Diehl教授合著重点讲解非线性最优控制问题的数值求解方法并配套基于CasADi与MPCTools的示例和习题便于在MATLAB环境中实操验证。资料为单本PDF电子书压缩包5.66MB共1个文件内容完整涵盖MPC理论基础、动态模型构建、优化问题设定、约束处理以及线性规划、二次规划、非线性规划等求解算法并包含参数整定、稳定性与鲁棒性分析等设计方法。已有376人学习下载适合作为MPC入门打基础与进阶研究的核心参考资料。1. 为什么第二版的 MPC 仍然值得你亲手写一遍从一个超调事故说起前年调试一台小型温控模拟项目X蒸汽阀和冷水阀互相耦合负载一切换温度就往上一冲PID 怎么压都压不住超调。后来把核心回路换成模型预测控制Model Predictive ControlMPC超调从 4℃ 压到 0.8℃约束、耦合和设定值斜坡一次性处理掉了。那是我第一次意识到MPC 不是只活在论文里的黑匣子而是一条能在现场复现的工程路径。“理论、计算、设计”三个词正好对应落地时的三道关怎么选预测模型、怎么把优化算出来、怎么接到真实装置上。这篇笔记不写书评只写我沿着这三条线把 MPC 跑通的方法参数怎么定、命令怎么写、约束怎么设、坑在哪里。做控制的人可以从这里直接上手。2. 模型选择与预测时域为什么设计决策比求解器更决定成败MPC 的起点不是求解器而是预测模型。书里大量公式落到工程现场最后往往只变成一个离散状态空间方程。这个方程的质量决定了后续所有调参工作有没有意义。本章先把模型选择说清楚再把预测时域、控制时域、采样周期和约束软硬这几个设计决策依次拆开。这些决策看起来是参数实际上是在定义优化问题的形状比换求解器影响大得多。2.1 预测模型怎么选线性状态空间是默认起点但不是唯一答案模型预测控制的核心不是“预测”本身而是用模型在优化器里推算未来。最常见也最容易上手的建模形式是线性状态空间import numpy as np Ts 0.1 A np.array([[1.0, Ts], [0.0, 1.0]]) B np.array([[0.0], [Ts]]) C np.array([[1.0, 0.0]])这是一个双积分器模型第一个状态是位置第二个状态是速度控制量是加速度。很多机械、液位和温控回路在小范围内都能近似成这个形式。A 阵里的1.0, Ts是前向欧拉离散的结果B 阵第二行是Ts表示控制量在一个采样周期内对速度的影响。C 阵只拿了位置做输出速度在真实系统里往往需要估计而不是直接测量。选择模型时我一般按这个顺序判断先看工作范围是否足够小如果是就用一个线性状态空间模型。如果工作范围大就用线性化得到多组模型做增益调度而不是一上来就上非线性 MPC。只有确实需要处理硬非线性比如化学反应速率、死区补偿才考虑非线性预测模型因为它每个采样周期要解一个非线性规划计算量和复杂度都高一个数量级。至于阶跃响应模型也就是动态矩阵控制那一类做法好处是不需要内部状态坏处是输出约束和观测器很难接。我只有在完全没有机理数据、只有现场阶跃响应曲线时才用它。这个选择不是理论洁癖而是决定后面解的是 QP 还是非线性规划、能不能上嵌入式、调参难度差多少。一句话线性状态空间是默认起点不要为了追求模型复杂度而堆公式。2.2 预测时域、控制时域和采样周期先定物理时间再定离散步数MPC 有几个基础数字每个都直接改变优化问题的大小和视野。我把常用的工程起点列在下面参数物理含义工程起点Ts控制周期主导时间常数的 1/10 到 1/20Np预测时域Np*Ts 覆盖主导时间常数的 3 到 5 倍Nc控制时域Np 的 1/5 到 1/2松弛变量权软约束惩罚强度1e3 到 1e5最常见的一个误区是直接定 Np而不看 Np 乘上 Ts 之后的物理时间。Np20、Ts0.1 意味着控制器往前看 2 秒如果 Ts 改成 0.2同样的 Np 就变成 4 秒。所以我的习惯是先算物理时间再反推步数。Ts 太小离散矩阵接近单位阵预测值对输入不敏感求解器数值条件变差Ts 太大一个扰动进来要等到下一个周期才响应。Np 也不是越大越好窗口太长前面几步的约束被后面大量无约束步数平均掉响应会变慢计算量按决策变量数量立方增长。Nc 越大控制器自由度越多但更容易让控制量高频振荡。我一般先把 Nc 放在 3 到 5 再调不要一开始就让它等于 Np。判断 Np 够不够有个很省事的办法把 Np 减半如果响应明显变差说明窗口不够把 Np 翻倍如果曲线几乎不变说明已经足够这时就回到原来的值。第二版教材里会花很大篇幅讨论终端约束和稳定性工程上如果 Np 覆盖了主导动态很多稳定性问题不会出现但如果 Np 太短加再多终端约束也难救。2.3 约束的硬软决策与代价函数权重先保可行再谈最优MPC 和 PID 最大的区别是能直接处理约束但约束怎么设决定它会不会“无解”。输入约束一定是硬约束因为执行器物理限位不能破。状态约束和输出约束建议用软约束预测有误差硬约束一旦冲突优化器返回 infeasible控制量直接断掉。软约束的标准做法是引入松弛变量 ε并给一个足够大的惩罚项eps cp.Variable(nonnegTrue) constraints.append(x[k 1] x_max eps) constraints.append(x[k 1] x_min - eps) cost 1e4 * eps**2这里的 ε 只惩罚预测越界不惩罚正常轨迹。权重 1e4 是经验起点太小会允许平时偷偷越界太大会退化成硬约束。代价函数本身我习惯写成状态误差的二次型加上控制量二次型。Q 决定哪些状态误差更“贵”R 决定控制量用得多贵。调权重时先看比例位置权重和速度权重的比值、状态权重和控制权重的比值而不是单个绝对值。如果位置是 0 到 100速度是 -1 到 1直接用 Q 的对角值去调速度项会被位置项淹没。先做归一化把所有状态除以各自的最大允许值再设 Q 和 R这样权重才真正反映优先级。还有一个容易踩的坑设定值跳变太大时硬约束会让优化问题在第一个周期就撞墙。我的习惯是给参考轨迹加一阶斜坡或低通滤波让设定值的变化时间略小于 Np*Ts。代价函数权重不是玄学但量纲和参考轨迹不处理好的话它看起来就是玄学。提示输入约束保持硬输出约束保持软。把这句话写成注释放在代码里能避免大量 infeasible 问题。3. 离散化与 QP 求解把 MPC 算出来的最小可运行代码MPC 从公式到可运行关键一步是把“预测”变成“优化”。所有状态和输入在预测时域里堆成向量代价函数变成二次型约束变成线性不等式于是每个采样周期要解一个二次规划。本章给出完整的最小实现并说明离散化方法和求解器参数怎么设。3.1 连续模型怎么离散零阶保持比前向欧拉更稳妥如果被控对象给的是连续模型第一步是离散化。常见做法是零阶保持认为输入在一个采样周期内保持不变这符合实际数模转换器的行为。用 Python 写很简单from scipy.signal import cont2discrete A_c np.array([[0.0, 1.0], [0.0, 0.0]]) B_c np.array([[0.0], [1.0]]) C_c np.array([[1.0, 0.0]]) D_c np.array([[0.0]]) Ts 0.1 sys_d cont2discrete((A_c, B_c, C_c, D_c), Ts, methodzoh) Ad, Bd, Cd, Dd sys_d[0], sys_d[1], sys_d[2], sys_d[3]这里cont2discrete用零阶保持做离散得到的 Ad 是矩阵指数Bd 是输入到状态的精确映射。前向欧拉只是用Ad I A*Ts在 Ts 足够小时误差不大但这里建议直接用 ZOH。注意B_c是 2 行 1 列因为系统只有一个输入D_c必须显式给出否则函数无法推断输出矩阵的维度。离散化完成后检查一下 Ad、Bd 的量级。如果 Ts 取得太大B 的数值会变大求解时容易出现数值问题。C 在 MPC 优化里不一定直接用但在后面的状态估计里一定会用到。3.2 滚动优化最小实现Python 加 cvxpy 搭一个通用 QP下面的代码是一个可以直接跑的线性 MPC。模型就用章节 2.1 的双积分器目标是把位置从 0 控制到 1速度归零。import numpy as np import cvxpy as cp Ts 0.1 A np.array([[1.0, Ts], [0.0, 1.0]]) B np.array([[0.0], [Ts]]) nx, nu B.shape Np 20 # 预测时域: 20 步, 对应 2 秒 Nc 5 # 控制时域: 前 5 步可以变化 Q np.diag([10.0, 0.1]) R np.array([[0.05]]) u_max, u_min 1.0, -1.0 x_max np.array([1.5, 2.0]) x_min -x_max def mpc_step(x0, ref): ref_vec np.array([ref, 0.0]) # 期望位置和期望速度 x cp.Variable((Np 1, nx)) u cp.Variable((Np, nu)) cost 0 constraints [x[0] x0] for k in range(Np): constraints.append(x[k 1] A x[k] B u[k]) constraints.append(u[k] u_max) constraints.append(u[k] u_min) constraints.append(x[k 1] x_max) constraints.append(x[k 1] x_min) cost cp.quad_form(x[k 1] - ref_vec, Q) cost cp.quad_form(u[k], R) # 控制时域之外的输入保持 Nc-1 步的值 for k in range(Nc, Np): constraints.append(u[k] u[Nc - 1]) prob cp.Problem(cp.Minimize(cost), constraints) prob.solve(solvercp.OSQP, verboseFalse, eps_abs1e-5, eps_rel1e-5, polishTrue) if prob.status ! cp.OPTIMAL: return None return u[0].value x np.array([0.0, 0.0]) for k in range(60): u_opt mpc_step(x, 1.0) if u_opt is None: print(infeasible at step, k) break x A x B np.array([u_opt[0]]) if k % 10 0: print(k, k, pos, x[0])这段代码把 MPC 的四个要素都放进去了。预测状态 x 和预测输入 u 是优化变量每个采样周期重新求解离散状态方程x[k1] A x[k] B u[k]被写成等式约束输入限幅是硬约束目标函数把跟踪误差和控制量放在同一个二次型里。Nc之后的输入被约束成最后一个可变值这是减少决策变量的常用做法等价于把未来输入参数化为分段常值。参数说明Np20 意味着预测窗口是 2 秒大于双积分模型从 0 到 1 的基本动态因此不需要额外终端约束也能稳定。Q 里位置权重 10、速度权重 0.1代表位置优先R0.05 表示控制量不太贵响应会比较快。想更保守把 R 调到 0.1 以上。polishTrue让 OSQP 在找到解后做一次抛光减少数值误差。如果求解状态不是 optimal函数返回 None控制器应做保持上一拍输出的保护逻辑而不是直接异常退出。3.3 求解器选型与热启动什么时候该把 QP 换成更大口径3.2 用的是 OSQP因为线性约束 MPC 的 QP 本质是凸优化。OSQP 是稀疏算子分裂方法对中小规模问题快且内存小几百个变量以内的模型用 CVXPY 封装很合适。如果预测模型带非线性就得换非线性求解器常见方案是 CasADi 配 IPOPT但每个周期求解时间不稳快回路要谨慎。如果目标是嵌入式或极小采样周期常见做法是把 QP 的矩阵结构和稀疏模式提前固化用代码生成工具导出 C 代码再在实时环境里调用这时就不在 Python 层解了。另一个实际技巧是热启动。MPC 相邻两个采样周期的问题几乎一样把上一时刻的最优解作为新问题的初值迭代次数能明显减少。在 CVXPY 里用warm_startTrue再把上一帧的u.value、x.value赋给变量。热启动的坑是问题变了但初值还不可行所以不要强行依赖它保留足够的max_iter和软约束。如果计算时间超过采样周期的一半就要警惕要么减小 Np 或 Nc要么把状态约束改成软约束要么把求解精度放宽到 1e-4。判断标准不是平均求解时间而是最坏情况下的求解时间MPC 必须保证每个周期都算完。4. 从仿真到实物接进 MPC 前先过模型失配、延迟和调参这三关代码能在仿真里跑离实物还差很远。我接实物前一定会做三件事故意让仿真对象和预测模型不一致把一拍延迟和输入变化率约束写进问题最后按固定顺序调参。这一章就是把这三关逐个说清楚。4.1 仿真怎么搭才有参考价值故意给对象模型掺一点假如果仿真对象用的 A/B 与 MPC 里的 A/B 完全一致那只是数学恒等式不是控制测试。真正有参考价值的仿真要主动制造模型失配。下面这段示意怎么搭rng np.random.default_rng(42) noise_std 0.01 def real_plant(x, u): # 对象模型: 增益故意比预测模型大 15% Ar np.array([[1.0, Ts], [0.0, 1.0]]) Br np.array([[0.0], [Ts * 1.15]]) return Ar x Br u x_true np.array([0.0, 0.0]) for k in range(100): u_opt mpc_step(x_est, 1.0) x_true real_plant(x_true, u_opt) y_meas x_true[0] rng.normal(0, noise_std) # x_est 在真实工程里由观测器更新real_plant的 Br 比预测模型的 B 大了 15%模拟阀门系数偏差和未建模摩擦。噪声加在测量上而不是状态上因为现场传感器的位置测量有噪声速度状态则需要估计。你会看到只调 Q/R 压不掉稳态误差因为控制器相信模型增益而模型增益是错的。真正的修正在下一节用观测器或增广扰动项。仿真要做的不是看起来好看而是主动制造失配去把问题压出来。仿真里不出现的坑实物上大概率加倍还回来。注意仿真对象和预测模型完全一致等于在证明同一个方程自己等于自己没有任何验证价值。4.2 一拍延迟、输入变化率与状态估计MPC 的隐藏前提MPC 的预测方程默认k 时刻测到状态k 时刻给定 ux[k1] 就更新。实物不是这样。实际链路是测量、计算、D/A、执行器响应至少要一拍时间才真正作用。常见做法有两种。第一种保持计算在采样周期中尽早完成u[k] 在 k 时刻直接下发只要计算时间小于 Ts 的 10%延迟可以忽略。第二种把延迟建模进预测增广一个状态u_d让u_d[k1] u[k]预测方程改成x[k1] A x[k] B u_d[k]。慢系统用第一种快系统或执行器响应慢的用第二种。输入变化率约束经常被漏掉。执行器不只有饱和还有速率限制。如果不限制相邻两步的变化量MPC 会在约束边界附近疯狂抖动。在约束里加du_max 0.2 constraints.append(u[0] - u_prev du_max) constraints.append(u[0] - u_prev -du_max) for k in range(1, Np): constraints.append(u[k] - u[k - 1] du_max) constraints.append(u[k] - u[k - 1] -du_max)u_prev是上一拍实际下发值。没有这一步OSQP 求出来的 u 可能是一拍 1、一拍 -1 的锯齿波R 小的时候尤其明显。状态估计更不用说MPC 要求全状态不是只有测量输出。位置可以直接测速度不能直接测时用 Kalman 滤波而不是简单差分。差分出来的速度噪声极大直接进 MPC 会让 u 抖这也是很多人第一次上 MPC 觉得控制器不稳定的主要原因。4.3 调参顺序与观察指标先看控制量再看跟踪曲线我一般按固定顺序调不要上来就动 Q/R。先把开环阶跃做了看主导时间常数和稳态增益然后按第 2 章的规则定 Ts、Np、Nc约束全部放开Q 和 R 按归一化量设Q 先只对角R 从中间值开始再把输入约束加回去观察求解时间和控制量抖动最后加状态约束、变化率约束和软约束做模型失配仿真。观察指标包括上升时间、超调量、控制量最大变化量、约束激活次数和最坏求解时间。指标期望超调量5% 以内控制量最大变化量不超过执行器速率上限的一半约束激活次数可预见、可接受最坏求解时间小于采样周期的 50%如果超调大先加大 Q 或减小 R如果抖动说明 R 太小或状态估计噪声大先别调 Q如果稳态误差不为零方向不是调参而是检查观测器和模型增益偏差。这些顺序我踩过不少坑最反直觉的是“超调大不一定是权重问题可能是 Np 太短”。调参前先把窗口拉长到主导时间常数的 5 倍再说很多问题自己就消失了。5. MPC 落地避坑五条让我返工最多的现场教训下面这些来自我调试过程中的记录每个现象都对应一个具体原因。有些看起来像玄学实际上是接口、时序或建模位置不对。按现象、原因、解决三条写清楚。5.1 现象一预测曲线漂亮实机一上去就抖现象模型仿真里曲线平滑控制量稳定一接到模拟执行器指令像噪声一样来回跳。原因最常见的是 R 太小或者 Nc 太大导致控制器把测量噪声和未建模动态当成需要快速补偿的偏差。另一个常见来源是没有状态估计直接用差分信号当速度。解决先看控制量曲线一个采样周期一跳基本就是上述问题。把 R 加大 5 到 10 倍Nc 减到 3再看曲线。如果加了低通滤波后控制量不再高频跳说明噪声确实从测量端进来。滤波时间常数不要超过 Ts/5否则会引入相位滞后控制器会反过来去补偿滞后。再不行就查执行器反馈和指令是否一致看有没有死区或回程间隙。5.2 现象二约束一激活就震荡甚至无解现象设定值接近边界时求解器返回 infeasible控制量保持上一拍值或者当 u 达到限幅后位置在边界附近来回冲。原因把输出或状态约束设成了硬约束。预测模型有任何误差预测路径就可能越过边界优化器找不到可行解。约束激活时震荡通常是参考设定值直接顶到边界控制器想在一个周期内追上去输入在上下限之间来回切换。解决输入约束保持硬输出约束改成软约束加松弛变量。给参考轨迹加斜坡让它知道边界不可瞬间到达。松弛变量权重取 1e4 左右只在快要越界时才起作用。出现 infeasible 时不要立刻放松输入约束先看是哪个约束在哪个预测步被违反把被违反的那一个改成软约束。5.3 现象三换一台设备就要重调模型参数看着全对现象同一套 MPC 代码和参数从一台设备搬到同型号另一台设备超调和调节时间差别明显。原因设备间增益和时间常数有分散性名义模型参数在开环辨识里看着合理但控制器对模型增益误差很敏感。增益偏差 10% 就足以让预测的稳态位置和目标错位闭环品质立刻变样。解决换设备后先做一次小阶跃把实测稳态增益用来修正 B 矩阵。现场有多个工况时把增益做成调度表而不是调 Q/R 去掩盖模型误差。增益调度后如果稳态误差仍在要在观测器里加扰动状态让 MPC 的预测模型拥有积分作用。5.4 现象四采样周期一换控制性能完全变样现象Ts 从 0.2 改成 0.1Np、Q/R 都没动响应速度明显变慢或者出现低频振荡。原因离散模型没有随 Ts 重新更新或者更新了但预测时域的物理时间变了。Np20 时Ts0.2 的窗口是 4 秒改成 0.1 后只有 2 秒控制器看得更短自然更保守。Nc5 也是一样物理上允许变化的时长从 1 秒变成 0.5 秒。解决每次改 Ts用cont2discrete重新生成 Ad/Bd。保持物理时间一致也就是同时调整 Np 和 Nc让 NpTs、NcTs 不变。改完之后重新看一眼开环阶跃确认预测模型的稳态增益没有因为离散化而变化。5.5 现象五代价函数权重调不动反而越调越差现象把 Q 里位置权重从 10 加到 1000闭环响应没有变快甚至出现求解精度下降和约束抖动。原因数量级差太多。位置是 0 到 100速度是 -1 到 1Q 里 1000 和 0.1 跨了 4 个数量级Hessian 条件数过高求解器在数值上分不清谁重要。另一个原因是参考值、状态都直接用工程单位权重被数值尺度吃掉了。解决先标幺化状态除以各自最大允许值输入除以执行器满量程。归一化后从 Qdiag([1, 0.1])、R[0.1] 开始调权重比例才有物理意义。调完后再换回工程单位但约束逻辑保持不变。这个坑最隐蔽因为它不报错只会让曲线越调越乱。5.6 现场排查顺序先查控制量再查状态估计最后查模型控制异常时我的排查顺序是固定的先打印控制量曲线和求解状态看是不是高频抖再看状态估计值是不是跟踪得上测量值最后才怀疑模型。求解状态出现非 optimal 时把prob.status打出来。如果总是 infeasible先把硬约束改成软约束把参考斜坡加长。状态估计和测量差得远先修观测器不要动 MPC 参数。模型失配导致的问题权重怎么调都没有用因为控制器在用错误的模型做预测。这个顺序能避免大多数越调越乱的情况。6. 进阶用闭环“两步检验法”验证你的 MPC 是否真正到位MPC 能在仿真里跑是一回事交付前我会固定做两道检验。第一步用与 MPC 完全一致的模型做闭环阶跃此时如果收敛慢或振荡说明基础窗口或权重有问题别急着甩锅给模型。第二步把对象模型的增益改 10% 到 20%加量测噪声再跑三个不同幅度的设定值观察超调、控制量抖动和无解次数。这一步过不了就回去改模型或约束而不是靠调参掩盖。具体判定指标超调小于 5%控制量最大变化量小于执行器速率上限的一半所有prob.status都是 optimal。如果第二步出现控制量抖优先检查 R 和状态估计如果稳态误差不为零优先检查扰动观测器和模型增益如果 infeasible优先检查软约束和参考斜坡。两步都通过后我会把参数放到一个独立配置块里Ts、Np、Nc、Q、R、约束、求解器容差全部集中管理换设备只改配置不碰控制代码。这个习惯帮我节省了大量翻找参数的时间也避免“上次在别的地方调过”这种说不清的状态。有一次我就是在双积分模型上把权重调得很漂亮忘记加模型失配和噪声结果一接实物就开始抖。随后补上两步检验法才把问题定位到状态估计而不是继续跟 Q/R 较劲。从那以后这一步成了我交付前的最后一道关卡。希望帮到你。本文还有配套的精品资源点击获取
