简介这是用于弹簧质量阻尼系统模型预测控制MPC仿真的MATLAB源码包适合学习控制理论、特别是二次型优化与预测控制算法的学生或工程师参考。压缩包共包含3个m文件分别承担系统建模、MPC控制器设计及单位阶跃响应仿真功能代码规模精简整体仅2KB便于逐行阅读和二次修改。作者已在MATLAB 2016a下实测运行附带的注释与参数设置能帮助读者理解从连续模型离散化、预测方程推导到标准二次型问题求解的完整流程并展示了解析法和数值法两种优化求解方式的实现差异。已有1028人学习下载作为小型教学示例具有较高参考价值适合与相关博文配合使用快速掌握MPC在机械系统中的落地写法。1. 弹簧质量阻尼系统的MPC控制仿真不调PID也能稳住的采样控制做运动控制的工程师多半被同一件事卡过弹簧质量阻尼系统这种二阶对象PID调得好好的一旦把采样周期拉长、或者非要硬约束阀门开度在0到1之间PID那套带宽和相位裕度的直觉就失灵了。弹簧质量阻尼系统SMD的MPC控制算法仿真m文件源码解决的就是这个矛盾——用一个显式的预测模型在每个采样周期里实时算出一段有限时域的最优控制动作把约束直接写进优化问题里。市面上大多数demo只给simulink框图真正能scilab式跑起来的m文件反而稀缺尤其是把状态空间离散化、QP求解器调用和约束矩阵拼装全部显式写出来的版本更是少见。这篇按我自己的落地套路来从连续模型转到离散状态空间再到MPC的代价函数与求解器选型最终给出可直接改参数运行的主程序和子函数最后聊参数整定和仿真发散的排错顺序。2. 先立模型弹簧质量阻尼系统的状态空间离散化与约束设定2.1 连续域微分方程到标准状态空间的换血过程弹簧质量阻尼系统SMD的物理方程是经典的二阶微分方程。设质量块质量为m弹簧刚度k阻尼系数c外力为u位移为x1x速度为x2ẋ那么牛顿第二定律给出m·x2_dot c·x2 k·x1 u大多数人写状态空间时会直接套标准形式但在做MPC之前需要把输入矩阵B的量纲理清楚。常见错误是把外力直接除以质量、然后塞进B矩阵导致后续预测序列全部差一个m倍。整理成标准的状态方程ẋ1 x2ẋ2 (u − k·x1 − c·x2) / m写成矩阵形式A_c [0, 1; −k/m, −c/m]B_c [0; 1/m]C_c [1, 0]这里是连续时间模型但数字控制器实际工作在离散域。连续A矩阵的指数映射是MPC预测的核心用零阶保持器ZOH离散化得到的离散状态空间才真正进入仿真代码。2.1.1 零阶保持器离散化的两种m文件写法第一种写法是直接用MATLAB的c2d函数这是最快路径几乎所有源码都会用。第二种是手写矩阵指数用expm(A_c·Ts)配合积分零阶保持公式A_d expm(A_c·Ts)B_d A_c(A_d − I)·B_c后者在只允许纯m文件、不允许依赖Control System Toolbox的环境里更稳妥。我一般会在源码开头放一个选择开关允许用c2d也可能用手写两个结果应该保持一致这本身就是一个验证点。延展一下MPC仿真里的离散周期Ts不是随便取的。弹簧质量阻尼系统的自然频率ωnsqrt(k/m)采样频率建议取10到20倍ωn再折算成Ts。比如k2m1ωn约1.414rad/sfs大约14到28Hz取Ts0.05s算是底线保守一点取0.02s。采样时间太短约束优化的计算量陡增太长则离散误差直接让预测模型失真。2.2 约束的物理含义与矩阵化装配代码MPC相较PID最大的价值是处理约束。在弹簧质量阻尼系统里约束主要来自执行器力u的上下限、变化率还有可能的质量块位移边界。源码里约束必须以线性不等式的形式进入QP求解器形式为lb ≤ u(k) ≤ ubu_min ≤ u(ki) − u(ki−1) ≤ u_max把控制序列向量化后约束条件就变成对U_all[u(k); u(k1); ...; u(kNc−1)]的取值空间做限制。翻译到m文件里就是拼装Ain矩阵和bin向量% 约束矩阵装配Nc为控制时域u_min/u_max为幅值上下限 Ain [eye(Nc); -eye(Nc)]; bin [u_max * ones(Nc,1); -u_min * ones(Nc,1)];这里Ain的维度是2Nc×Nc上半块约束上界下半块约束下界。注意负号的位置非常容易被搞混。更严格的做法是把状态约束也叠进去。输出位移x1的硬约束格式是x1_min ≤ C_d·x(ki) ≤ x1_max但状态约束在MPC里属于是“软硬兼施”的问题硬约束存在可行解为空的可能专门做约束软化slack variable会在后续版本升级中体现。考虑变化率约束时再拼一行对角差分矩阵复杂度可控但不是所有源码都覆盖。2.2.1 模型失配时的约束可行性考量弹簧质量阻尼系统的k和c参数在大部分教材题目中是恒定的但真实系统里弹簧可能磨损、阻尼系数随温度漂移。MPC的“闭环鲁棒性”依赖于模型预测准确模型失配会让约束边界失效。做仿真时给参数加±10%的摄动再跑一遍如果约束被持续违反说明Np和Nc的搭配需要重新考虑。3. MPC核心算法流程代价函数构造与m文件求解器调用逻辑3.1 预测模型在m文件里的迭代形式有了离散状态空间A_d/B_d/C_d之后MPC算法的预测部分就是在当前状态x(k)的基础上把未来Np步的状态x(k1)到x(kNp)推演出来。注意这里Np是预测时域Nc是控制时域通常Nc≤Np超过Nc之后的控制量保持u(kNc−1)不变。预测方程的矩阵化是MPC仿真里提升速度的关键也是源码里的必考知识点。把所有预测状态堆成一个大列向量X_all把所有控制量堆成U_all预测可以一步写成X_allF_mat·x(k)G_mat·U_all。其中F_mat和G_mat由A_d与B_d反复迭代生成两个矩阵的维度分别是Np·nx×(nx)和Np·nx×Nc。用手写双重循环当然可行但矩阵化更快且代码更短更重要的是和推导公式一一对应。规格是nx代表状态维度对弹簧质量阻尼系统即2nu代表输入维度1Np长度预测时域Nc控制时域。% 预测矩阵构造状态维度nx2输入维度nu1 F_mat zeros(Np * nx, nx); G_mat zeros(Np * nx, Nc); for i 1:Np F_mat((i-1)*nx1:i*nx, :) A_d^i; for j 1:min(i, Nc) G_mat((i-1)*nx1:i*nx, (j-1)*nu1:j*nu) A_d^(i-j) * B_d; end end逻辑说明F_mat里的A_d^i表示从当前状态到第i步的开环状态转移G_mat每个块A_d^(i−j)·B_d表示第j步的控制输入对第i步预测状态的贡献只有j≤i时有意义j≥i的块保持零。这正好对应因果性控制量只能影响未来时刻的状态。若用NcNp循环内min(i,Nc)保证了超出控制时域的部分没有新增控制输入参与。3.2 代价函数展开与Q、R权重矩的阵设定标准MPC代价函数写成J Σ(x(ki)^T·Q·x(ki)) Σ(u(kj)^T·R·u(kj))端值惩罚项terminal cost在调节问题里也常加入但大部分基础仿真源码会省略末尾阶段这种做法在短预测时域时稳定裕度会缩小在无限时域最优控制理论里端值项对应着离散黎卡提方程的解想严谨可以做扩展。把x的预测表达式代入代价函数后J就变成U_all的二次函数J U_all^T·H·U_all 2·x(k)^T·f^T·U_all constant% 代价函数矩阵H和线性项fQ_state为状态权重R_in为输入权重 H G_mat * kron(eye(Np), Q_state) * G_mat kron(eye(Nc), R_in); f_mat G_mat * kron(eye(Np), Q_state) * F_mat * x_now;H是Nc×Nc的对称正定矩阵Q正半定、R正定时成立f_mat是Nc×1。这个二次规划的标准形式quadprog要求f是列向量维度对不上就会报错。kron函数把Q_state的作用扩展到整个预测时域每个时刻的状态误差都计入代价这是MPC和LQR一脉相承的地方。3.2.1 Q和R的物理指向弹簧质量阻尼系统里Q阵的对角元分别对应位移误差和速度误差的惩罚。Q(1,1)调大则位移收敛更激进代价是控制量更大Q(2,2)调大则速度衰减更猛容易造成超调。R对应控制力的大小惩罚R增大系统变得“懒散”。初始仿真建议Qdiag([10,1])或Qdiag([100,10])R1Np20Nc3这样一组中庸参数起步。3.3 m文件中调用quadprog还是自己写求解器在MATLAB仿真里quadprog是标配。它支持’interior-point-convex’、’active-set’等算法约束多时active-set收敛更快约束少或没有约束时直接用闭式解更高效。对弹簧质量阻尼系统的仿真无约束情况可以直接用解析解U*−H(f_mat)可以绕过求解器大幅提升仿真速度。有约束的情况下调用方式如下options optimoptions(quadprog, Display, off, Algorithm, active-set); [U_opt, ~, exitflag] quadprog(H, f_mat, Aineq, bineq, [], [], [], [], [], options); u_applied U_opt(1);逻辑说明U_opt是Nc×1的优化结果按照MPC滚动优化原则只取第一个分量u_applied作为当前实际控制量下一时刻重新测量状态再求解一轮。exitflag需要检查正数表示成功收敛负数可能是无解或迭代超限。约束包含上下界时quadprog的lb/ub参数直接给不需要重复拼装进Aineq但提前用Aineq形式写一遍便于控制时域内不同段取不同界。3.3.1 无约束解析解与有约束数值解的切换建议仿真调试初期先跑无约束版本验证预测模型的正确性再打开约束。这样分离了“模型错误”和“约束冲突”两类故障。无约束调试时用U*−H\f_mat注意H必须正定否则求解结果震荡甚至发散这是矩阵病态的一个典型预兆。调试顺序建议模型开环验证 → 无约束MPC → 只有幅值约束 → 幅值变化率约束4. 弹簧质量阻尼系统仿真m文件落地主程序、子函数与参数调优4.1 主程序结构与m文件目录组织一个可运行的MPC仿真m文件至少拆成三块主脚本、MPC控制器函数、被控对象模型函数。把模型和控制器分开写在两个子函数里后续换对象模型比如换成双质量弹簧系统不用动控制器。主脚本的流程是初始化参数 → 生成离散模型 → 构造预测矩阵 → 循环仿真每步求解QP再更新状态 → 画图。弹簧质量阻尼系统仿真m文件源码的核心循环如下% 主仿真循环总时长20s步长Ts0.05s Nsim 400; x_hist zeros(2, Nsim1); u_hist zeros(1, Nsim); x_now x0; for k 1:Nsim F_mat build_F(Np, A_d, nx); G_mat build_G(Np, Nc, A_d, B_d, nx, nu); U_opt solve_mpc(H, f_mat, Aineq, bineq, x_now); u_applied U_opt(1); u_hist(k) u_applied; x_now A_d * x_now B_d * u_applied; x_hist(:, k1) x_now; end在这个循环里F_mat和G_mat如果在每次迭代重新构造会拖慢仿真速度标准做法是在循环外一次生成前提是A_d和B_d不变化。x_now是当前状态的真实值在纯仿真中用上一时刻的A_d·xB_d·u递推在硬件在环时必须替换成传感器测量值。4.1.1 build_F和build_G子函数的完整m文件代码function F_mat build_F(Np, A_d, nx) F_mat zeros(Np * nx, nx); for i 1:Np F_mat((i-1)*nx1:i*nx, :) A_d^i; end end function G_mat build_G(Np, Nc, A_d, B_d, nx, nu) G_mat zeros(Np * nx, Nc * nu); for i 1:Np for j 1:min(i, Nc) G_mat((i-1)*nx1:i*nx, (j-1)*nu1:j*nu) A_d^(i-j) * B_d; end end end参数说明nx是状态数2nu是控制数1Np和Nc是时域参数。A_d^i用矩阵幂运算当i0时MATLAB返回单位阵正好对应d步前控制量的直接传递。这个函数设计的核心是让i和j的循环边界保持一致性只要错一个索引预测就会发散的莫名其妙。4.2 输出响应曲线与稳定性评估方法仿真结束后看三张图x1位置时间响应、u控制力时间响应、状态相平面图。位置响应曲线应该平滑收敛到目标值默认调节问题目标x0或阶跃目标控制力在约束边界内。数据层面计算收敛时间和超调量比如位置响应超调超过20%说明Q权重偏低或Np不够长。还有一个评估点是稳态误差。MPC在没有积分作用的情况下对模型失配会产生稳态偏差。从这个角度考虑弹簧质量阻尼系统仿真多数是调节问题目标是原点模型失配时不加积分项会有小残差。如果要把目标从0改到xref非零值状态方程写法换成x̃x−xref控制量也做平移变换不要把目标直接硬塞给代价函数。4.2.1 阶跃响应的目标平移操作示例% 目标位移xref误差状态xe x - [xref; 0] xe_now x_now - [xref; 0]; [U_opt, ~, flag] quadprog(H, f_mat * xe_now, Aineq, bineq);逻辑说明代价函数里的线性项f_mat乘上误差状态xe_now意味着预测模型追踪目标xref的偏移量控制量u_applied不需要额外补偿因为系统模型包含输入的直接作用项只要误差收敛到零控制量自然到达对应稳态值。此方案在无模型失配时零稳态误差。4.3 参数整定顺序Np、Nc、Q/R到底先动哪个MPC调试新手喜欢一上来就动Q/R矩阵重量其实效果最小。正确的调参顺序首先是Np预测时域其次是Nc控制时域最后才调权重。Np太小比如2、3时预测信息太少控制效果接近高增益比例控制容易出现震荡Np太大50以上时远期的预测几乎没有信息量还徒增计算耗时。弹簧质量阻尼系统这类振荡环节Np至少要覆盖一个振荡周期以上的时长按照周期T≈2π/ωd来计算。以k2m1c0.2为例系统阻尼比低于临界振荡周期大约4.4s取Ts0.05s时Np应该≥90才稳妥。但实际仿真里Np30就已经能工作得很好原因是最优控制本身不要求预测覆盖整个衰减过程只要求当前窗口内有合理的梯度方向。Nc的调法更简单从1起步逐步增大到某个阈值后几乎不影响响应曲线这时停止。经验起始参数Np50, Nc5, Qdiag([10;1]), R15. 仿真发散的排查次序与一个自动寻参技巧5.1 发散来源分类模型错误、离散化错误和QP无解仿真发散是MPC落地最常见的坑而且现象五花八门有的位置量直接NaN有的控制量振荡幅度越来越大有的仿真到一半报quadprog无解。排查顺序有讲究。第一排查对象是A_d和B_d是否组装正确。做法是单独跑一段开环x(k1)在不加控制的情况下应该是衰减或震荡趋势如果一仿真就指数型发散且不受约束控制那A_d的符号一定错了。测试命令很简单给一个初位移用离散模型递推100步不施加控制看是否稳定。第二排查H矩阵是否正定。无约束MPC用H\f_mat求解析解时如果H近奇异求解结果的幅值会异常大。正定性检验看eig(H)是否全为正注意Q_state为零矩阵时H只是半正定数值解仍然有可能工作但理论保证弱化。第三排查约束相互矛盾导致QP无解。比如此时位移超出硬约束边界同时控制力又在极限处即使把控制量打到最大也无法把位移拉回约束内。quadprog的exitflag输出为负值时先看约束边界是否合理再考虑加约束软化。5.1.1 收敛状态下消耗时间的量化验证弹簧质量阻尼系统的MPC仿真性能瓶颈绝大多数在quadprog的迭代求解。无约束解析解大概几毫秒有约束QP可能要几十毫秒到数百毫秒不等。用tic/toc包裹主循环统计总耗时如果显式Q的维度Nc·nu为10以上、耗时超过1秒先怀疑约束数量是否过载再考虑是否把内部iteration限制。5.2 基于闭环仿真的Np自动扫描m文件调参数靠手工一个个试效率太低且Np和Nc之间常有耦合效应。写一个简单的循环扫描脚本在Np数组和目标收敛时间之间建立映射自动选出最优参数。这也是很多源码附带的价值点% 扫描Np从10到100记录位置响应收敛到5%误差范围内的时间 Np_list 10:10:100; settle_time zeros(size(Np_list)); for idx 1:length(Np_list) Np Np_list(idx); % 重新构建矩阵H和f_mat并跑一次仿真 [~, settle_time(idx)] run_smd_mpc(Np, Nc, Q_state, R_in); endrun_smd_mpc把主循环封装成函数返回位置响应序列和稳定时间。这个扫描在几分钟内能完成远远快于手工试参。权重Q/R的扫描稍复杂因为Q和R的相对比例才是关键简单做法是把R固定为1Q改成q_param·eye(2)再扫描q_param。同时需要在代价函数里加入控制变化量的惩罚以压低执行器磨损公式扩展成J Σ x^T·Q·x Σ u^T·R·u Σ Δu^T·S·ΔuΔuu(k)−u(k−1)进入代价后H矩阵多出与时域差分矩阵相关的一项结构仍然保持二次型。当前时刻Δu中包含了上一次的控制量u(k−1)需要在f_mat中额外补一项与历史状态相关的偏差向量。这样弹簧质量阻尼系统的MPC仿真就不只停留在“能跑”的层面可以直接过渡到执行器有磨损限制的工程场景。本文还有配套的精品资源点击获取
弹簧质量阻尼系统MPC控制仿真：从状态空间离散化到QP求解的m文件实现
简介这是用于弹簧质量阻尼系统模型预测控制MPC仿真的MATLAB源码包适合学习控制理论、特别是二次型优化与预测控制算法的学生或工程师参考。压缩包共包含3个m文件分别承担系统建模、MPC控制器设计及单位阶跃响应仿真功能代码规模精简整体仅2KB便于逐行阅读和二次修改。作者已在MATLAB 2016a下实测运行附带的注释与参数设置能帮助读者理解从连续模型离散化、预测方程推导到标准二次型问题求解的完整流程并展示了解析法和数值法两种优化求解方式的实现差异。已有1028人学习下载作为小型教学示例具有较高参考价值适合与相关博文配合使用快速掌握MPC在机械系统中的落地写法。1. 弹簧质量阻尼系统的MPC控制仿真不调PID也能稳住的采样控制做运动控制的工程师多半被同一件事卡过弹簧质量阻尼系统这种二阶对象PID调得好好的一旦把采样周期拉长、或者非要硬约束阀门开度在0到1之间PID那套带宽和相位裕度的直觉就失灵了。弹簧质量阻尼系统SMD的MPC控制算法仿真m文件源码解决的就是这个矛盾——用一个显式的预测模型在每个采样周期里实时算出一段有限时域的最优控制动作把约束直接写进优化问题里。市面上大多数demo只给simulink框图真正能scilab式跑起来的m文件反而稀缺尤其是把状态空间离散化、QP求解器调用和约束矩阵拼装全部显式写出来的版本更是少见。这篇按我自己的落地套路来从连续模型转到离散状态空间再到MPC的代价函数与求解器选型最终给出可直接改参数运行的主程序和子函数最后聊参数整定和仿真发散的排错顺序。2. 先立模型弹簧质量阻尼系统的状态空间离散化与约束设定2.1 连续域微分方程到标准状态空间的换血过程弹簧质量阻尼系统SMD的物理方程是经典的二阶微分方程。设质量块质量为m弹簧刚度k阻尼系数c外力为u位移为x1x速度为x2ẋ那么牛顿第二定律给出m·x2_dot c·x2 k·x1 u大多数人写状态空间时会直接套标准形式但在做MPC之前需要把输入矩阵B的量纲理清楚。常见错误是把外力直接除以质量、然后塞进B矩阵导致后续预测序列全部差一个m倍。整理成标准的状态方程ẋ1 x2ẋ2 (u − k·x1 − c·x2) / m写成矩阵形式A_c [0, 1; −k/m, −c/m]B_c [0; 1/m]C_c [1, 0]这里是连续时间模型但数字控制器实际工作在离散域。连续A矩阵的指数映射是MPC预测的核心用零阶保持器ZOH离散化得到的离散状态空间才真正进入仿真代码。2.1.1 零阶保持器离散化的两种m文件写法第一种写法是直接用MATLAB的c2d函数这是最快路径几乎所有源码都会用。第二种是手写矩阵指数用expm(A_c·Ts)配合积分零阶保持公式A_d expm(A_c·Ts)B_d A_c(A_d − I)·B_c后者在只允许纯m文件、不允许依赖Control System Toolbox的环境里更稳妥。我一般会在源码开头放一个选择开关允许用c2d也可能用手写两个结果应该保持一致这本身就是一个验证点。延展一下MPC仿真里的离散周期Ts不是随便取的。弹簧质量阻尼系统的自然频率ωnsqrt(k/m)采样频率建议取10到20倍ωn再折算成Ts。比如k2m1ωn约1.414rad/sfs大约14到28Hz取Ts0.05s算是底线保守一点取0.02s。采样时间太短约束优化的计算量陡增太长则离散误差直接让预测模型失真。2.2 约束的物理含义与矩阵化装配代码MPC相较PID最大的价值是处理约束。在弹簧质量阻尼系统里约束主要来自执行器力u的上下限、变化率还有可能的质量块位移边界。源码里约束必须以线性不等式的形式进入QP求解器形式为lb ≤ u(k) ≤ ubu_min ≤ u(ki) − u(ki−1) ≤ u_max把控制序列向量化后约束条件就变成对U_all[u(k); u(k1); ...; u(kNc−1)]的取值空间做限制。翻译到m文件里就是拼装Ain矩阵和bin向量% 约束矩阵装配Nc为控制时域u_min/u_max为幅值上下限 Ain [eye(Nc); -eye(Nc)]; bin [u_max * ones(Nc,1); -u_min * ones(Nc,1)];这里Ain的维度是2Nc×Nc上半块约束上界下半块约束下界。注意负号的位置非常容易被搞混。更严格的做法是把状态约束也叠进去。输出位移x1的硬约束格式是x1_min ≤ C_d·x(ki) ≤ x1_max但状态约束在MPC里属于是“软硬兼施”的问题硬约束存在可行解为空的可能专门做约束软化slack variable会在后续版本升级中体现。考虑变化率约束时再拼一行对角差分矩阵复杂度可控但不是所有源码都覆盖。2.2.1 模型失配时的约束可行性考量弹簧质量阻尼系统的k和c参数在大部分教材题目中是恒定的但真实系统里弹簧可能磨损、阻尼系数随温度漂移。MPC的“闭环鲁棒性”依赖于模型预测准确模型失配会让约束边界失效。做仿真时给参数加±10%的摄动再跑一遍如果约束被持续违反说明Np和Nc的搭配需要重新考虑。3. MPC核心算法流程代价函数构造与m文件求解器调用逻辑3.1 预测模型在m文件里的迭代形式有了离散状态空间A_d/B_d/C_d之后MPC算法的预测部分就是在当前状态x(k)的基础上把未来Np步的状态x(k1)到x(kNp)推演出来。注意这里Np是预测时域Nc是控制时域通常Nc≤Np超过Nc之后的控制量保持u(kNc−1)不变。预测方程的矩阵化是MPC仿真里提升速度的关键也是源码里的必考知识点。把所有预测状态堆成一个大列向量X_all把所有控制量堆成U_all预测可以一步写成X_allF_mat·x(k)G_mat·U_all。其中F_mat和G_mat由A_d与B_d反复迭代生成两个矩阵的维度分别是Np·nx×(nx)和Np·nx×Nc。用手写双重循环当然可行但矩阵化更快且代码更短更重要的是和推导公式一一对应。规格是nx代表状态维度对弹簧质量阻尼系统即2nu代表输入维度1Np长度预测时域Nc控制时域。% 预测矩阵构造状态维度nx2输入维度nu1 F_mat zeros(Np * nx, nx); G_mat zeros(Np * nx, Nc); for i 1:Np F_mat((i-1)*nx1:i*nx, :) A_d^i; for j 1:min(i, Nc) G_mat((i-1)*nx1:i*nx, (j-1)*nu1:j*nu) A_d^(i-j) * B_d; end end逻辑说明F_mat里的A_d^i表示从当前状态到第i步的开环状态转移G_mat每个块A_d^(i−j)·B_d表示第j步的控制输入对第i步预测状态的贡献只有j≤i时有意义j≥i的块保持零。这正好对应因果性控制量只能影响未来时刻的状态。若用NcNp循环内min(i,Nc)保证了超出控制时域的部分没有新增控制输入参与。3.2 代价函数展开与Q、R权重矩的阵设定标准MPC代价函数写成J Σ(x(ki)^T·Q·x(ki)) Σ(u(kj)^T·R·u(kj))端值惩罚项terminal cost在调节问题里也常加入但大部分基础仿真源码会省略末尾阶段这种做法在短预测时域时稳定裕度会缩小在无限时域最优控制理论里端值项对应着离散黎卡提方程的解想严谨可以做扩展。把x的预测表达式代入代价函数后J就变成U_all的二次函数J U_all^T·H·U_all 2·x(k)^T·f^T·U_all constant% 代价函数矩阵H和线性项fQ_state为状态权重R_in为输入权重 H G_mat * kron(eye(Np), Q_state) * G_mat kron(eye(Nc), R_in); f_mat G_mat * kron(eye(Np), Q_state) * F_mat * x_now;H是Nc×Nc的对称正定矩阵Q正半定、R正定时成立f_mat是Nc×1。这个二次规划的标准形式quadprog要求f是列向量维度对不上就会报错。kron函数把Q_state的作用扩展到整个预测时域每个时刻的状态误差都计入代价这是MPC和LQR一脉相承的地方。3.2.1 Q和R的物理指向弹簧质量阻尼系统里Q阵的对角元分别对应位移误差和速度误差的惩罚。Q(1,1)调大则位移收敛更激进代价是控制量更大Q(2,2)调大则速度衰减更猛容易造成超调。R对应控制力的大小惩罚R增大系统变得“懒散”。初始仿真建议Qdiag([10,1])或Qdiag([100,10])R1Np20Nc3这样一组中庸参数起步。3.3 m文件中调用quadprog还是自己写求解器在MATLAB仿真里quadprog是标配。它支持’interior-point-convex’、’active-set’等算法约束多时active-set收敛更快约束少或没有约束时直接用闭式解更高效。对弹簧质量阻尼系统的仿真无约束情况可以直接用解析解U*−H(f_mat)可以绕过求解器大幅提升仿真速度。有约束的情况下调用方式如下options optimoptions(quadprog, Display, off, Algorithm, active-set); [U_opt, ~, exitflag] quadprog(H, f_mat, Aineq, bineq, [], [], [], [], [], options); u_applied U_opt(1);逻辑说明U_opt是Nc×1的优化结果按照MPC滚动优化原则只取第一个分量u_applied作为当前实际控制量下一时刻重新测量状态再求解一轮。exitflag需要检查正数表示成功收敛负数可能是无解或迭代超限。约束包含上下界时quadprog的lb/ub参数直接给不需要重复拼装进Aineq但提前用Aineq形式写一遍便于控制时域内不同段取不同界。3.3.1 无约束解析解与有约束数值解的切换建议仿真调试初期先跑无约束版本验证预测模型的正确性再打开约束。这样分离了“模型错误”和“约束冲突”两类故障。无约束调试时用U*−H\f_mat注意H必须正定否则求解结果震荡甚至发散这是矩阵病态的一个典型预兆。调试顺序建议模型开环验证 → 无约束MPC → 只有幅值约束 → 幅值变化率约束4. 弹簧质量阻尼系统仿真m文件落地主程序、子函数与参数调优4.1 主程序结构与m文件目录组织一个可运行的MPC仿真m文件至少拆成三块主脚本、MPC控制器函数、被控对象模型函数。把模型和控制器分开写在两个子函数里后续换对象模型比如换成双质量弹簧系统不用动控制器。主脚本的流程是初始化参数 → 生成离散模型 → 构造预测矩阵 → 循环仿真每步求解QP再更新状态 → 画图。弹簧质量阻尼系统仿真m文件源码的核心循环如下% 主仿真循环总时长20s步长Ts0.05s Nsim 400; x_hist zeros(2, Nsim1); u_hist zeros(1, Nsim); x_now x0; for k 1:Nsim F_mat build_F(Np, A_d, nx); G_mat build_G(Np, Nc, A_d, B_d, nx, nu); U_opt solve_mpc(H, f_mat, Aineq, bineq, x_now); u_applied U_opt(1); u_hist(k) u_applied; x_now A_d * x_now B_d * u_applied; x_hist(:, k1) x_now; end在这个循环里F_mat和G_mat如果在每次迭代重新构造会拖慢仿真速度标准做法是在循环外一次生成前提是A_d和B_d不变化。x_now是当前状态的真实值在纯仿真中用上一时刻的A_d·xB_d·u递推在硬件在环时必须替换成传感器测量值。4.1.1 build_F和build_G子函数的完整m文件代码function F_mat build_F(Np, A_d, nx) F_mat zeros(Np * nx, nx); for i 1:Np F_mat((i-1)*nx1:i*nx, :) A_d^i; end end function G_mat build_G(Np, Nc, A_d, B_d, nx, nu) G_mat zeros(Np * nx, Nc * nu); for i 1:Np for j 1:min(i, Nc) G_mat((i-1)*nx1:i*nx, (j-1)*nu1:j*nu) A_d^(i-j) * B_d; end end end参数说明nx是状态数2nu是控制数1Np和Nc是时域参数。A_d^i用矩阵幂运算当i0时MATLAB返回单位阵正好对应d步前控制量的直接传递。这个函数设计的核心是让i和j的循环边界保持一致性只要错一个索引预测就会发散的莫名其妙。4.2 输出响应曲线与稳定性评估方法仿真结束后看三张图x1位置时间响应、u控制力时间响应、状态相平面图。位置响应曲线应该平滑收敛到目标值默认调节问题目标x0或阶跃目标控制力在约束边界内。数据层面计算收敛时间和超调量比如位置响应超调超过20%说明Q权重偏低或Np不够长。还有一个评估点是稳态误差。MPC在没有积分作用的情况下对模型失配会产生稳态偏差。从这个角度考虑弹簧质量阻尼系统仿真多数是调节问题目标是原点模型失配时不加积分项会有小残差。如果要把目标从0改到xref非零值状态方程写法换成x̃x−xref控制量也做平移变换不要把目标直接硬塞给代价函数。4.2.1 阶跃响应的目标平移操作示例% 目标位移xref误差状态xe x - [xref; 0] xe_now x_now - [xref; 0]; [U_opt, ~, flag] quadprog(H, f_mat * xe_now, Aineq, bineq);逻辑说明代价函数里的线性项f_mat乘上误差状态xe_now意味着预测模型追踪目标xref的偏移量控制量u_applied不需要额外补偿因为系统模型包含输入的直接作用项只要误差收敛到零控制量自然到达对应稳态值。此方案在无模型失配时零稳态误差。4.3 参数整定顺序Np、Nc、Q/R到底先动哪个MPC调试新手喜欢一上来就动Q/R矩阵重量其实效果最小。正确的调参顺序首先是Np预测时域其次是Nc控制时域最后才调权重。Np太小比如2、3时预测信息太少控制效果接近高增益比例控制容易出现震荡Np太大50以上时远期的预测几乎没有信息量还徒增计算耗时。弹簧质量阻尼系统这类振荡环节Np至少要覆盖一个振荡周期以上的时长按照周期T≈2π/ωd来计算。以k2m1c0.2为例系统阻尼比低于临界振荡周期大约4.4s取Ts0.05s时Np应该≥90才稳妥。但实际仿真里Np30就已经能工作得很好原因是最优控制本身不要求预测覆盖整个衰减过程只要求当前窗口内有合理的梯度方向。Nc的调法更简单从1起步逐步增大到某个阈值后几乎不影响响应曲线这时停止。经验起始参数Np50, Nc5, Qdiag([10;1]), R15. 仿真发散的排查次序与一个自动寻参技巧5.1 发散来源分类模型错误、离散化错误和QP无解仿真发散是MPC落地最常见的坑而且现象五花八门有的位置量直接NaN有的控制量振荡幅度越来越大有的仿真到一半报quadprog无解。排查顺序有讲究。第一排查对象是A_d和B_d是否组装正确。做法是单独跑一段开环x(k1)在不加控制的情况下应该是衰减或震荡趋势如果一仿真就指数型发散且不受约束控制那A_d的符号一定错了。测试命令很简单给一个初位移用离散模型递推100步不施加控制看是否稳定。第二排查H矩阵是否正定。无约束MPC用H\f_mat求解析解时如果H近奇异求解结果的幅值会异常大。正定性检验看eig(H)是否全为正注意Q_state为零矩阵时H只是半正定数值解仍然有可能工作但理论保证弱化。第三排查约束相互矛盾导致QP无解。比如此时位移超出硬约束边界同时控制力又在极限处即使把控制量打到最大也无法把位移拉回约束内。quadprog的exitflag输出为负值时先看约束边界是否合理再考虑加约束软化。5.1.1 收敛状态下消耗时间的量化验证弹簧质量阻尼系统的MPC仿真性能瓶颈绝大多数在quadprog的迭代求解。无约束解析解大概几毫秒有约束QP可能要几十毫秒到数百毫秒不等。用tic/toc包裹主循环统计总耗时如果显式Q的维度Nc·nu为10以上、耗时超过1秒先怀疑约束数量是否过载再考虑是否把内部iteration限制。5.2 基于闭环仿真的Np自动扫描m文件调参数靠手工一个个试效率太低且Np和Nc之间常有耦合效应。写一个简单的循环扫描脚本在Np数组和目标收敛时间之间建立映射自动选出最优参数。这也是很多源码附带的价值点% 扫描Np从10到100记录位置响应收敛到5%误差范围内的时间 Np_list 10:10:100; settle_time zeros(size(Np_list)); for idx 1:length(Np_list) Np Np_list(idx); % 重新构建矩阵H和f_mat并跑一次仿真 [~, settle_time(idx)] run_smd_mpc(Np, Nc, Q_state, R_in); endrun_smd_mpc把主循环封装成函数返回位置响应序列和稳定时间。这个扫描在几分钟内能完成远远快于手工试参。权重Q/R的扫描稍复杂因为Q和R的相对比例才是关键简单做法是把R固定为1Q改成q_param·eye(2)再扫描q_param。同时需要在代价函数里加入控制变化量的惩罚以压低执行器磨损公式扩展成J Σ x^T·Q·x Σ u^T·R·u Σ Δu^T·S·ΔuΔuu(k)−u(k−1)进入代价后H矩阵多出与时域差分矩阵相关的一项结构仍然保持二次型。当前时刻Δu中包含了上一次的控制量u(k−1)需要在f_mat中额外补一项与历史状态相关的偏差向量。这样弹簧质量阻尼系统的MPC仿真就不只停留在“能跑”的层面可以直接过渡到执行器有磨损限制的工程场景。本文还有配套的精品资源点击获取