简介基于CasADi与IPOPT求解非线性模型预测控制NMPC的完整Matlab实现聚焦车辆/发动机轨迹跟踪与动力学控制场景适合自动化、计算机、电子信息、数学等专业的学生用于课程设计、期末大作业或毕业设计。代码采用参数化编程核心参数可灵活修改注释详细思路清晰并附有多种可直接运行的案例数据可快速体验从建模、求解到结果可视化的完整流程。包内共21个文件以15个m脚本为主覆盖主控程序、动力学模型、Pacejka轮胎魔术公式、RBF/辅助函数以及动画测试等模块另有csv、xlsx真实发动机数据、mat保存的仿真结果和txt说明文档便于对照运行与二次开发。压缩包整体约841KB轻量易用。目前已有414人学习适合希望深入掌握CasADi/IPOPT接口调用和NMPC工程实现的进阶学习者。 做非线性模型预测控制这几年我发现一个挺尴尬的现状理论大家都懂代价函数怎么选、约束怎么加也都清楚但真正想让NMPC在MATLAB里变成一个能仿真、能调试、能往实物上搬的东西卡住的往往是工具链。网上能找到的代码要么是线性MPC要么是教学版的单输入单输出真到了强非线性、强耦合的系统线性MPC压不住自写的梯度求解器又不稳定。这套基于CasADi和IPOPT的MATLAB NMPC代码正好把最难的那一环补上了——用符号框架把非线性预测控制问题完整建模再用内点法求解器稳定地解出来。这篇文章我把自己读代码、跑通、调参的过程完整记录下来适合正在做NMPC仿真、准备做嵌入式部署或者刚接触CasADi工具链的同学参考。1. 为什么这套代码值得细看NMPC的工程定位1.1 线性MPC够用为什么还要非线性很多工控场景MPC的线性版本就已经能干活了。水箱液位、温度调节、简单的速度控制模型在工作点附近线性化后的误差可以接受。但一旦系统跑在大范围工况切换里或者动力学本身带强非线性例如摆杆大幅度摆动、机械臂变构型、反应器温度跨区间线性模型的预测就会在几个步长之后严重跑偏。MPC的命根子就是“预测”模型不准后面优化的结果再好也是错的。非线性模型预测控制在每个控制周期不做事先的线性化而是直接用非线性微分方程去预测未来一段时间的状态轨迹然后在非线性约束下求解一个最优控制问题。代价是计算量变大好处是模型始终反映真实系统的行为。尤其处理强耦合、输入受限、状态受限的场合NMPC的优势是碾压级的。这套代码最大的价值就是用CasADi的自动微分和符号优化能力把非线性预测的建模和求解门槛降了下来。1.2 CasADi和IPOPT在整个链路中的角色第一次看到CasADi这个名字容易把它当成又一个优化求解器。实际上它不是一个求解器而是一个符号数学与优化建模框架。你可以用CasADi定义符号变量、构造非线性函数、做自动微分然后把这些表达式直接喂给底层的数值求解器。NMPC里最麻烦的雅可比矩阵和海森矩阵计算CasADi能自动推导而且是精确的解析导数不是差分近似这对IPOPT收敛的帮助巨大。IPOPT则是真正干苦力活的内点法非线性求解器。它接受CasADi吐出来的目标函数、约束条件和导数信息在每次迭代中沿着障碍函数路径找最优点。对一般规模的NMPC问题控制在几十个状态变量、几十步预测时域之内IPOPT的收敛速度和稳健性都属于第一梯队。再加上MATLAB版CasADi有mex接口数据传递效率不差整个链路在桌面端做实时仿真完全可行。2. 搭建NMPC的整体思路与建模2.1 用SX符号变量描述系统模型在写任何代价函数和约束之前首先要把被控对象的数学模型变成CasADi认识的符号表达式。这套代码用的是SX类型也就是用稀疏符号矩阵来承载状态变量、输入变量和系统动态。之所以不直接用普通的MATLAB函数是因为优化过程中需要反复求梯度符号表达式的自动微分能力在这个环节能省一大笔事。import casadi.* % 状态x1是角度x2是角速度 x SX.sym(x, 2); u SX.sym(u, 1); % 被控对象参数单摆 m 0.2; % 质量 L 0.3; % 摆长 b 0.02; % 阻尼系数 g 9.81; % 连续时间动力学x_dot f(x, u) x_dot [x(2); -(g/L)*sin(x(1)) - (b/(m*L^2))*x(2) u/(m*L^2)]; f Function(f, {x, u}, {x_dot});这里用单摆做例子主要是为了方便讲清楚逻辑实际工程系统换成机械臂、四旋翼、反应器也一样只要能把状态方程写出来流程完全一致。写到这一步CasADi已经知道系统的结构和每个导数怎么算了。2.2 RK4离散化从连续方程到预测模型NMPC的预测模型必须能向前推演未来轨迹连续微分方程不能直接用于离散的优化计算所以要先做数值离散化。这套代码用的是经典四阶龙格库塔也就是RK4。RK4单步精度高对大多数机械系统来说控制周期在10到50毫秒时离散误差可以完全忽略。Ts 0.05; % 控制周期 % 单步RK4积分把连续模型f变成离散映射F(x,u) - x_next Xk x; for k 1:4 k1 f(Xk, u); k2 f(Xk Ts/2*k1, u); k3 f(Xk Ts/2*k2, u); k4 f(Xk Ts*k3, u); Xk Xk Ts/6*(k1 2*k2 2*k3 k4); end F Function(F, {x, u}, {Xk});注意这里的F是一个新的CasADi函数对象输入是当前状态和当前输入输出是下一时刻状态。这个离散映射就是NMPC闭环仿真里的“预测步”在优化问题的每个步长上都会反复调用它。用RK4替代简单的欧拉法虽然计算多了一点但可以让我们用更大的预测步长而保持预测精度这点在后端求解效率上的收益很划算。3. 核心代码逐段解读单摆跟踪控制3.1 定义优化变量和参量预测模型准备好了接下来就是构造离散时间最优控制问题。这段代码是这个项目里最需要细看的部分。N 20; % 预测时域 opti casadi.Opti(); % 创建优化问题容器 % 决策变量状态轨迹和控制轨迹 X opti.variable(2, N1); U opti.variable(1, N); % 参数变量当前测量状态每次迭代更新 P opti.parameter(2, 1); % 参考轨迹这里用定值目标 x_ref [pi/3; 0];opti.variable定义出来的矩阵就是这个NLP的决策变量包含整个预测时域内的所有状态和所有控制输入。用Opti栈写NMPC的好处是思路直接你再也不用手动把变量整理成一维向量、维护索引关系了CasADi内部帮你全部打理好。这里P是parameter它的含义是“外部传入的参数”在优化问题构造时数值未知求解前再set_value进去。每个控制周期我们把当前系统的测量状态赋给P然后重新求解这一点是理解NMPC滚动优化的关键。3.2 代价函数与约束构造代价函数决定控制器“在乎什么”。这个例子里用目标跟踪误差的二次型和输入能量的二次型也是工程上最常用的套路。Q diag([10; 0.1]); % 状态误差权重 R 0.01; % 输入权重 cost 0; for k 1:N err X(:,k) - x_ref; cost cost err * Q * err U(k) * R * U(k); end opti.minimize(cost);Q矩阵里第一个元素是角度误差权重给得大控制器会优先把角度拉回目标值第二个元素是角速度权重给得小允许一定的动态过程。R控制输入的激进程度太小会让控制量猛冲太大则响应变慢。这些权重在模型正确的前提下直接决定了闭环动态品质。约束方面一是系统动力学约束把每个步长之间的状态演化关系写死opti.subject_to(X(:,1) P); % 初始状态等于当前测量 for k 1:N opti.subject_to(X(:,k1) F(X(:,k), U(:,k))); end二是物理约束比如输入转矩不能超过极限、状态不能越过安全范围opti.subject_to(-5 U 5); % 输入饱和 opti.subject_to(-pi X(1,:) pi); % 角度限幅输入饱和约束非常关键。如果模型假设控制量没有上限优化算法会利用这个漏洞输出一个匪夷所思的大信号导致真实系统里执行器根本跟不上去闭环性能瞬间崩塌。3.3 求解配置和闭环仿真问题定义好接下来声明求解器类型。这里用的是IPOPT求解阶段把所有信息交给它。opti.solver(ipopt, struct(print_time, 0), struct(print_level, 0));print_time和print_level都设为0是为了让求解过程安静一点避免每个控制周期都刷屏在实时仿真里这很有用。接下来是闭环仿真主循环% 系统初始状态 x0 [0; 0]; x_current x0; sim_steps 200; for t 1:sim_steps % 传入当前测量状态 opti.set_value(P, x_current); % 初始猜测用上一步的解做warm start opti.set_initial(X, repmat(x_ref, 1, N1)); opti.set_initial(U, zeros(1, N)); % 求解 sol opti.solve(); % 提取控制序列并施加第一个控制量 u_applied sol.value(U(1)); % 用真实系统模型推进这里用高精度积分模拟被控对象 [~, y_next] ode45((t, y) pendulum_dynamics(y, u_applied), [0 Ts], x_current); x_current y_next(end, :); % 记录状态和控制轨迹用于画图 history.x(:, t) x_current; history.u(t) u_applied; end这里有一个细节在真实工程中系统下一步的状态来自传感器测量或状态估计器而不是我们代码里自己积分的模型。仿真里用ode45推进是为了模拟一个“真实对象”而控制器内部使用的是RK4离散模型两者存在误差这正好检验了NMPC对模型误差的鲁棒性。3.4 warm start从慢到快的钥匙如果每次求解都从零开始IPOPT需要较多迭代才能找到可行解。NMPC里一个成熟技巧是用上一时刻的最优解作为当前时刻的初始猜测也就是warm start。上一控制周期的状态轨迹和控制序列与当前周期应当非常接近从那里出发能让求解器在前几次迭代就落到一个次优但可行的区域迭代次数大幅下降。% 上一个求解周期的解 prev_X sol.value(X); prev_U sol.value(U); % 移一位作为新周期的初始猜测 opti.set_initial(X, [prev_X(:,2:end), prev_X(:,end)]); opti.set_initial(U, [prev_U(:,2:end), prev_U(:,end)]);这个细节我认为是这套代码里最有工程价值的部分之一。刚开始用CasADi做NMPC时我每次都忽略初始猜测导致求解时间忽长忽短实时性完全没法保证。加上warm start之后平均每步求解时间下降了将近一半而且时间波动小了很多。4. IPOPT求解器配置与调参经验4.1 求解器参数怎么调才稳IPOPT是一套参数很多的内点法求解器但对NMPC场景常调的其实就那么几个。我用下来最关键的是线搜索容差tol、可接受的终止容差acceptable_tol以及最大迭代次数max_iter。options struct(... print_level, 0, ... max_iter, 200, ... tol, 1e-4, ... acceptable_tol, 1e-6, ... acceptable_iter, 5, ... mu_strategy, adaptive, ... linear_solver, mumps); opti.solver(ipopt, struct(), options);tol设置的是标准收敛阈值acceptable_tol是放宽版的收敛标准。有些工况下NMPC并不需要把最优解算到小数点后六位控制目标达到1e-4级别就足够这时候acceptable_tol就能让我们提前终止迭代把计算时间让给下一个控制周期。mu_strategy建议设成adaptive。默认的内点法屏障参数更新策略有时候会让求解器在可行域边界来回试探自适应策略收敛更平滑。线性求解器我默认用mumps因为它随IPOPT一起发布不需要额外配置。如果系统规模大或者遇到数值困难可以考虑ma57但需要商业许可证。4.2 时域长度、权重和约束的工程取舍预测时域N的选择是NMPC调试里最需要反复试的一环。这个例子里N20对应1秒的预测长度。对小系统来说这个长度一般够用但如果被控对象动态响应慢比如温度控制预测时域需要覆盖系统大部分瞬态N要相应加大。代价是决策变量变多求解时间非线性上涨。我常用的判断方法是给系统一个阶跃输入看状态到达稳态的调整时间预测时域大约是调整时间的1到1.5倍。如果N太长而计算资源紧张优先减少控制输入序列的自由度比如让控制量在最后几步保持恒定这是一个在精度和实时性之间找平衡的常用折中方案。权重矩阵方面Q和R的比值比绝对值更重要。Q中状态误差权重和R中输入权重的比值决定了控制系统的“刚柔”程度。调试时先固定R逐步放大Q观察闭环响应达到可接受速度后停止。不要在第一步就追求完美参数NMPC的强项在于模型预测过大的Q容易导致约束边界处剧烈抖动这是一种典型的过调表现。4.3 实时性分析瓶颈在哪里在MATLAB里跑NMPC求解器计算只占一部分时间。每个控制周期里将当前状态设置进参数、处理变量映射、提取解这些接口调用也有固定开销。实测下来15到30步预测时域、2到4维状态的小规模问题在普通PC上单步求解时间大约在20到100毫秒量级这个速度对秒级或毫秒级的控制系统来说已经能用。如果要做高速实时控制建议往两个方向走一是减少N或改用更宽松的容差但这会牺牲控制品质二是用CasADi的代码生成功能把整个优化问题编译成C代码脱离MATLAB的解析开销。代码生成是CasADi最厉害的特性之一同一套NLP定义可以原封不动地发给CodeGenerator生成一个独立的C函数再编译成mex或共享库。我在某个实际项目里用这个方法把求解时间压缩了大约4倍效果非常显著。5. 常见问题与调试实录5.1 求解失败“Infeasible Problem Detected”这是NMPC新手最容易撞上的报错。IPOPT在问题无可行解时会直接给出类似提示。出现这个报错九成原因是约束设置自相矛盾。最容易忽略的是状态限幅和动力学约束之间的关系比如预测时域太短系统物理上不可能在几步之内从一个状态跑到另一个满足限幅的状态优化器就只能罢工。排查思路很简单先把所有不等式约束松弛掉只保留动力学等式约束运行看是否可行。如果能收敛说明问题出在约束边界上逐步把约束加回来每一步都运行测试定位到具体是哪个约束导致不可行。工程上还有一个更稳妥的方案是引入松弛变量把硬的状态约束变成软约束这在CasADi里的实现非常方便。5.2 闭环发散但求解器报告“成功”比求解失败更烦人的是IPOPT每次都成功收敛但闭环仿真结果飘了。这种情况大部分是离散化精度和预测模型与仿真模型不一致造成的。我遇到过类似的问题控制器内部用RK4离散积分而仿真端用欧拉法模拟被控对象控制周期取大了之后预测轨迹和真实轨迹迅速分道扬镳导致控制器陷入“盲人骑瞎马”的窘境。另一个常见原因是状态约束过紧导致的极限环。系统为了满足约束在每个周期内来回调整控制量但更新率跟不上系统动态就会出现发散。解决办法是检查约束是否合理特别是状态限幅是否小于系统的自然工作范围。如果约束本身没问题试试减小时间步长这会同时改善预测精度和求解稳定性。5.3 MATLAB环境下的CasADi细节在MATLAB里用CasADi有个容易踩坑的地方CasADi的Opti对象在循环内反复调用solve()默认情况下解的状态不会自动保留需要手动通过set_initial重新赋值。如果没有每步更新初始猜测IPOPT会从零开始求解时间会明显增加。这个我在3.3小节里强调过实际操作时一定别省这一步。另一个细节是符号变量和数值类型转换。CasADi的SX变量在求解器返回后是DM类型不能直接用于MATLAB矩阵运算要用full()或者value()转成普通数值。别小看这个转换嵌套在复杂仿真里时类型不匹配的报错会让人找半天。写代码时建议在闭环接口处就做一次彻底的转换不要在后续运算里再用DM类型做加减乘除。还有一点是版本兼容性。CasADi的MATLAB接口会跟某个具体MATLAB版本绑定编译下载时一定要认准和你的MATLAB版本匹配的发布包。曾经有同行在R2024a环境里装了旧版CasADi结果mex文件直接加载失败浪费了大半天排查环境问题。装好之后第一件事跑一下官方自带的test函数确认接口正常再开始写代码。6. 一些运行代码时的实操细节解压这套代码后我建议按“先跑通、再改模型、后调参数”的顺序推进。先别着急改代码把原始的示例直接运行一遍确认你的MATLAB版本和CasADi版本兼容、mex文件可以正常加载、IPOPT求解器能求解成功。很多人在第一步环境上就卡住了后面再漂亮的模型也是纸上谈兵。跑通之后再做三类改动。第一类是把模型换成本身的仿真对象这是NMPC代码迁移的核心凡是涉及系统方程的位置都要检查包括连续模型定义、离散化函数以及仿真端的被控对象模型。第二类是修改代价函数和约束根据实际控制目标调整Q、R矩阵和上下限。第三类是调性能参数在求解时间与预测精度之间找到适合自己的平衡点。另外一个小技巧是给闭环仿真记录观测数据把状态轨迹、控制输入、每一步的求解时间和IPOPT的退出状态存下来分析。退出状态尤为重要IPOPT返回的终止标志会明确告诉你是因为收敛条件满足退出还是因为最大迭代次数用尽退出。后者意味着求解器可能返回了一个并非最优的解但NMPC的滚动优化结构对这种次优解有一定的容忍度只要实时性优先可以适当放宽收敛标准前提是你要知道自己放弃了什么。这套代码的另一个可以扩展的方向是加入状态估计环节。现在它假设状态完全可测这在实际系统里往往不成立。可以在NMPC前面接一个扩展卡尔曼滤波器或者Moving Horizon Estimator而MHE本身也是CasADi的强项同一套符号框架可以同时处理估计和控制想想就知道这是多么顺手的组合。我个人的经验是NMPC这东西第一遍跑通是难点但只要工具链用熟了后续的迭代和优化速度会非常快而CasADi加上IPOPT正是这条路上最值得投入学习的一套方案。本文还有配套的精品资源点击获取