1. 项目概述从“黑箱”到“白箱”的控制艺术在工业自动化、机器人、自动驾驶这些领域我们总希望系统能“看得远一点”别等到撞墙了才刹车。这就是模型预测控制Model Predictive Control, MPC的核心魅力——它不像传统的PID那样只盯着当前的误差“头痛医头”MPC更像一个老练的棋手能基于对未来的几步推演走出当前最优的一步。这个项目就是带你亲手用Matlab这把“瑞士军刀”把MPC从理论公式变成一行行可运行的代码去驾驭离散、连续、线性乃至非线性的各种“模型马匹”。简单说MPC的工作流程是一个滚动优化的闭环在每个控制周期它利用一个描述系统动态的数学模型预测未来一段时间内系统状态的变化然后求解一个优化问题找到一序列最优的控制输入使得预测轨迹最贴合我们的期望比如快速、平稳、省能量最后只实施这序列中的第一个控制量。到了下一个周期重复这个过程用最新的测量值刷新预测这就是“滚动时域”的精髓。所以MPC的性能高度依赖于两个东西预测模型的准确度以及优化问题的求解效率。我们这个项目要做的就是针对不同类型的模型离散/连续、线性/非线性在Matlab里搭建起这个“预测-优化”的框架。为什么用Matlab因为它提供了一个从模型建立、控制器设计、仿真验证到性能分析的完整生态。无论是用ss、tf定义线性状态空间用ode45解算非线性微分方程还是用quadprog、fmincon求解优化问题Matlab都有成熟的工具箱和清晰的语法能让我们把精力集中在控制逻辑本身而不是底层算法的实现上。接下来我们就从最基础的线性离散模型开始一步步深入到更复杂的场景。2. 核心思路与方案选型如何为你的系统挑选“模型水晶球”在动手写代码之前我们必须想清楚我的系统到底该用哪种模型来描述这直接决定了后续控制器设计的复杂度和计算负担。这里没有银弹只有权衡。2.1 离散 vs. 连续时间观的抉择这个选择关乎你如何看待时间。连续时间模型通常由微分方程描述比如dx/dt Ax Bu。它认为时间是连续流动的最能反映物理世界的本质。在Matlab里我们常用Simulink或者用ode系列函数来仿真这类系统。离散时间模型通常由差分方程描述比如x(k1) A_d * x(k) B_d * u(k)。它认为时间是被采样成一个个瞬间的这与数字控制器的执行方式每隔一个采样周期Ts计算一次天然契合。关键决策点对于MPC设计我们几乎总是使用离散时间模型。原因很直接MPC的优化是在每个离散的采样时刻进行的预测的未来状态也是基于离散时间步的。即便你的物理系统是连续的也需要先将其离散化。在Matlab中对于线性系统你可以用c2d函数轻松将连续模型(A, B, C, D)转换为离散模型(A_d, B_d, C_d, D_d)并指定采样周期Ts。对于非线性系统离散化可能需要用到欧拉法、龙格-库塔法等数值方法。2.2 线性 vs. 非线性复杂度的分水岭这是MPC实现难度最大的分界线。线性MPC (LMPC)系统模型是线性的优化问题通常是二次型目标函数线性约束是凸的可以转化为**二次规划QP**问题。Matlab的Model Predictive Control Toolbox就是为此而生或者你也可以用quadprog手动构建。它的求解速度极快适合高速实时控制。非线性MPC (NMPC)系统模型是非线性的优化问题通常是非凸的需要求解**非线性规划NLP**问题。这要复杂得多计算量巨大常用fmincon这类求解器。虽然强大但实时性是个挑战。实操心得不要一上来就追求NMPC。首先尝试用线性模型。很多非线性系统可以在工作点附近进行线性化用线性模型设计MPC往往能取得不错的效果。如果非线性很强或者工作范围很广再考虑NMPC。在Matlab中对于连续非线性模型dx/dt f(x, u)可以在平衡点(x0, u0)处使用jacobian或linmod函数进行线性化。2.3 方案选型矩阵为了更直观我们可以根据系统特性快速定位系统特性推荐模型类型Matlab核心工具/函数优点缺点/挑战动态平滑工作点固定线性离散模型ss,c2d,mpc,quadprog设计简单求解快理论成熟对非线性、大范围动态性能不佳动态平滑工作点变化多模型/增益调度线性MPC多个mpc对象在线切换比单一线性模型适应范围广切换逻辑需精心设计动态非线性强但可求导非线性连续模型离散化后用于NMPCode45(仿真),fmincon(优化)模型精度高控制性能潜力大计算量大实时性差可能陷入局部最优模型复杂或为黑箱基于数据的线性模型如ARX或简单非线性模型arx,nlarx,System Identification Toolbox无需精确机理模型依赖数据质量外推能力弱我们这个项目将覆盖前三种情况给出具体的Matlab实现范例。3. 线性离散MPC的Matlab手把手实现我们从最经典、应用最广的线性离散MPC开始。假设我们已经有了一个离散线性状态空间模型x(k1) A * x(k) B * u(k)y(k) C * x(k)我们的目标是让输出y跟踪参考信号r同时控制输入u不要变化太剧烈。3.1 问题构建把控制目标写成数学公式MPC在每个时刻k要解的优化问题通常长这样Minimize J Σ [ (y(ki|k) - r(ki))^T * Q * (y(ki|k) - r(ki)) ] Σ [ Δu(ki|k)^T * R * Δu(ki|k) ] i1..Hp i0..Hc-1 Subject to: x(ki1|k) A * x(ki|k) B * u(ki|k) y(ki|k) C * x(ki|k) u_min ≤ u(ki|k) ≤ u_max Δu_min ≤ Δu(ki|k) ≤ Δu_maxHp: 预测时域看多远。Hc: 控制时域未来多少步控制量可以优化通常Hc ≤ HpHc之后控制量保持不变。Q,R: 权重矩阵Q调跟踪性能R调控制力度。Q越大跟踪越紧R越大控制动作越柔和。Δu(k) u(k) - u(k-1)是对控制量增量的约束能使系统更平稳。3.2 代码实现从零搭建QP求解框架我们不直接调用mpc工具箱而是用quadprog手动构建这能让你彻底理解MPC的“内脏”。%% 1. 定义系统模型 (离散双积分器示例位置控制) Ts 0.1; % 采样时间 A [1 Ts; 0 1]; % 状态矩阵 [位置; 速度] B [0.5*Ts^2; Ts]; % 输入矩阵 C [1 0]; % 输出矩阵输出位置 D 0; sysd ss(A, B, C, D, Ts); nx size(A,1); % 状态维度 nu size(B,2); % 输入维度 ny size(C,1); % 输出维度 %% 2. 设置MPC参数 Hp 20; % 预测时域 Hc 5; % 控制时域 Q 10 * eye(ny); % 输出误差权重 R 0.1 * eye(nu); % 控制增量权重 % 约束 u_min -2; u_max 2; du_min -1; du_max 1; %% 3. 构建QP问题的矩阵 (这是核心步骤) % 3.1 构建预测矩阵 [Phi, Gamma] build_prediction_matrices(A, B, C, Hp, Hc); % 需要自定义这个函数 % 3.2 构建QP标准形式: min (1/2) * z * H * z f * z, s.t. Aineq * z bineq % 决策变量 z [Δu(k|k); Δu(k1|k); ...; Δu(kHc-1|k)] H 2 * (Gamma * kron(eye(Hp), Q) * Gamma kron(eye(Hc), R)); % f 与参考轨迹和当前状态有关在每一步循环中计算 % 3.3 构建输入和输入增量的约束矩阵 [Aineq_u, bineq_u] build_input_constraints(Hc, nu, u_min, u_max, du_min, du_max); % 需要自定义 % 可能还有状态/输出约束... %% 4. 模拟闭环控制 Nsteps 100; x zeros(nx, Nsteps1); u zeros(nu, Nsteps); y zeros(ny, Nsteps); r 1.0 * ones(ny, Nsteps); % 阶跃参考信号 x_k [0; 0]; % 初始状态 u_km1 0; % 上一时刻控制量 for k 1:Nsteps % 计算QP中的 f 向量 % F (Y_ref - Ψ * x_k) * Qbar * Γ其中Qbar kron(eye(Hp), Q) Y_ref repmat(r(:,k), Hp, 1); % 假设未来参考值不变 Psi Phi(1:ny*Hp, :); % 输出预测中与初始状态相关的部分 f -2 * (Y_ref - Psi * x_k) * kron(eye(Hp), Q) * Gamma; % 构建当前时刻的输入约束 bineq % 需要将 u_min, u_max 转换为关于 Δu 的约束并考虑 u(k-1) bineq_u_k bineq_u; bineq_u_k(1:nu) bineq_u_k(1:nu) u_km1; % 修正上界约束中的常数项 bineq_u_k(nu1:2*nu) bineq_u_k(nu1:2*nu) - u_km1; % 修正下界约束中的常数项 % 求解QP options optimoptions(quadprog, Display, off); [z_opt, ~, exitflag] quadprog(H, f, Aineq_u, bineq_u_k, [], [], [], [], [], options); if exitflag 0 du_k z_opt(1:nu); % 取第一个控制增量 u_k u_km1 du_k; else warning(QP求解失败保持上一时刻控制量); u_k u_km1; end % 施加控制更新系统状态 (这里用真实模型模拟) x_k A * x_k B * u_k; y_k C * x_k; % 存储数据 x(:,k1) x_k; u(:,k) u_k; y(:,k) y_k; % 为下一步更新 u_km1 u_k; end %% 5. 绘图 figure; subplot(2,1,1); plot(1:Nsteps, y, b-, 1:Nsteps, r(1,1:Nsteps), r--); legend(输出y, 参考r); title(输出跟踪性能); subplot(2,1,2); stairs(1:Nsteps, u); title(控制输入u);上面代码中的build_prediction_matrices和build_input_constraints是两个关键的自定义函数它们负责将MPC预测方程转化为QP标准形式。这是线性MPC实现中最需要细心推导的部分。注意事项手动构建QP矩阵时矩阵维度非常容易出错。一个有效的调试方法是先用一个非常简单的系统比如一阶系统和小时域Hp3, Hc2来验证你的矩阵构建函数。打印出Phi、Gamma矩阵手动推算一两个预测步看看是否吻合。4. 非线性MPC的Matlab实现探索当系统模型为x(k1) f(x(k), u(k))离散非线性或dx/dt f(x, u)连续非线性时我们就进入了NMPC的领域。这里我们展示基于连续非线性模型、通过直接转录法离散化后使用fmincon求解的思路。4.1 直接转录法把连续问题“打散”成离散问题我们不再像线性MPC那样推导全局的预测矩阵而是将整个预测时域上的状态和控制量都作为优化变量并让动力学方程作为约束条件。%% 非线性MPC示例控制一个简单倒立摆到竖直位置 % 模型连续非线性状态为 [角度θ; 角速度θ_dot]控制为力矩u % 动力学方程 d/dt [θ; θ_dot] [θ_dot; (g/l)*sin(θ) u/(m*l^2)] %% 1. 定义参数和模型函数 m 1; l 1; g 9.81; odefun (t, x, u) [x(2); (g/l)*sin(x(1)) u/(m*l*l)]; %% 2. 设置NMPC参数 Ts 0.05; % 控制/离散化周期 Hp 25; % 预测时域 N Hp; % 离散点数 nx 2; nu 1; % 权重 Q diag([10, 1]); % 状态误差权重 R 0.01; % 控制量权重 % 约束 u_min -5; u_max 5; %% 3. 定义NLP的目标函数和约束函数 % 决策变量 Z [x0; u0; x1; u1; ...; x_{N-1}; u_{N-1}; x_N] % 其中 xi 是状态ui 是控制量。x0 是当前状态参数也是优化的起点。 total_vars N*(nxnu) nx; % 变量总数 % 目标函数最小化预测时域内的状态偏差和控制能量 cost_func (Z) nmpc_cost(Z, N, nx, nu, Q, R); % 非线性约束动力学方程通过欧拉离散化 nonlcon (Z) nmpc_constraints(Z, N, nx, nu, Ts, odefun); % 边界约束 lb -inf(total_vars,1); ub inf(total_vars,1); % 控制量约束 for i 0:N-1 idx_u 1 i*(nxnu) nx; % ui的位置 lb(idx_u) u_min; ub(idx_u) u_max; end %% 4. 模拟闭环 sim_time 5; time_steps floor(sim_time / Ts); x_history zeros(nx, time_steps1); u_history zeros(nu, time_steps); x0 [0.2; 0]; % 初始状态稍偏离零点 Z_initial_guess zeros(total_vars, 1); % 初始猜测可以更智能 options optimoptions(fmincon, Algorithm,interior-point, ... MaxIterations, 100, Display, off, ... SpecifyObjectiveGradient,false, ... SpecifyConstraintGradient,false); % 用数值梯度 for k 1:time_steps % 将当前状态 x0 作为参数固定到决策变量前部 Z_initial_guess(1:nx) x0; % 求解NLP [Z_opt, ~, exitflag] fmincon(cost_func, Z_initial_guess, ... [], [], [], [], lb, ub, nonlcon, options); if exitflag 0 warning(fmincon未收敛使用上一时刻解或零输入); u_k 0; else % 提取第一个控制量 u_k Z_opt(1 nx); % Z [x0; u0; ...] end % 模拟系统到下一时刻 (使用更精确的积分器如ode45) [~, x_traj] ode45((t,x) odefun(t,x,u_k), [0 Ts], x0); x0 x_traj(end, :); % 存储 x_history(:, k1) x0; u_history(:, k) u_k; % 为下一次优化准备初始猜测采用“热启动”将本次解平移作为下次初猜 Z_initial_guess warm_start_shift(Z_opt, N, nx, nu); end这里省略了nmpc_cost,nmpc_constraints和warm_start_shift这三个函数的详细代码。它们分别负责计算目标函数遍历离散点累加(xi - x_ref)*Q*(xi - x_ref) ui*R*ui。施加动力学约束对于每个离散区间要求x_{i1} x_i Ts * f(x_i, u_i)欧拉法计算约束违反量。热启动将本次最优解[x0, u0, x1, u1, ..., xN]平移用[x1, u1, x2, u2, ..., xN, u_{N-1}]作为下一次优化的初始猜测能极大加速收敛。实操心得与避坑指南计算负担NMPC的决策变量规模是(Hp * (状态数控制量数))fmincon求解可能很慢不适合高动态或采样周期很短的系统。工业上常用实时迭代RTI或显式NMPC等高级方法。离散化方法欧拉法最简单但精度差可能导致仿真不稳定。对于刚性或精度要求高的系统应考虑在约束函数中使用更精确的离散化方法如梯形法、龙格-库塔法或者采用正交配点法等更专业的直接转录法。初始猜测至关重要一个糟糕的初始猜测如全零可能导致fmincon收敛到局部最优甚至不收敛。“热启动”是提升NMPC在线计算效率最有效的技巧之一。梯度提供如果能为fmincon提供目标函数和约束的解析梯度甚至Hessian矩阵求解速度和可靠性会大幅提升。可以使用符号计算工具箱或自动微分来生成这些导数。5. 工具箱快速上手Matlab MPC Toolbox对于工业应用和快速原型Matlab自带的Model Predictive Control Toolbox是更稳健和便捷的选择。它完美支持线性时不变LTI模型的MPC设计。%% 使用MPC工具箱设计控制器 % 假设已有离散线性模型 sysd (来自前面示例) % 1. 创建MPC对象 Ts_mpc Ts; % 采样时间需与模型一致 p 20; % 预测时域 m 5; % 控制时域 mpcobj mpc(sysd, Ts_mpc, p, m); % 2. 配置控制器参数 % 设置权重 mpcobj.Weights.OutputVariables [10]; % 输出权重 (Q) mpcobj.Weights.ManipulatedVariablesRate 0.1; % 控制增量权重 (R) % 设置约束 mpcobj.MV.Min u_min; mpcobj.MV.Max u_max; mpcobj.MV.RateMin du_min; mpcobj.MV.RateMax du_max; % 设置输出信号类型跟踪参考 mpcobj.Model.Nominal.Y 0; % 可设置工作点 % 3. 进行闭环仿真 Tf 10; % 仿真总时间 r 1; % 参考信号 simopt mpcsimopt(); simopt.RefLookAhead off; % 参考预览关闭 simopt.Constraints on; simopt.OpenLoop off; % 使用 sim 命令仿真 [~, ~, ~, info] sim(mpcobj, Tf/Ts, r, [], simopt); % info 结构体包含了详细的优化信息 % 4. 更常见的用法在Simulink中使用MPC Controller模块 % 将mpcobj导出到工作空间在Simulink中拖入MPC Controller模块 % 在模块参数中指定控制器对象为mpcobj并连接好参考信号、输出反馈和控制输出。注意事项MPC Toolbox功能强大但“黑箱”程度较高。在将控制器投入实际应用前务必使用mpcmove命令或Simulink进行充分的闭环仿真测试特别是测试在模型失配、测量噪声和外部扰动下的鲁棒性。工具箱还提供了sensitivity、review(mpcobj)等函数用于分析控制器性能。6. 常见问题、调试技巧与性能优化在实际实现和调试MPC时你会遇到各种各样的问题。下面是一些典型问题及其排查思路。6.1 QP求解失败或无可行解问题quadprog返回exitflag为负值或无解。排查检查约束是否过紧这是最常见原因。特别是控制增量约束du_min/du_max和控制量约束u_min/u_max可能存在冲突。例如系统当前需要一个大控制量来纠正误差但du_max限制它只能缓慢增加导致在预测时域内永远无法达到所需控制量。临时放宽约束或调整权重是诊断方法。检查预测模型是否正确确保你构建的Phi和Gamma矩阵是正确的。用一个已知的初始状态和控制序列手动计算几步预测与矩阵乘法结果对比。检查QP矩阵的正定性矩阵H必须是半正定的。确保你的Q和R权重矩阵是半正定和正定的。6.2 控制器性能不佳震荡、超调大、响应慢问题输出跟踪有较大超调或持续震荡。调参指南调整权重Q和R这是最直接的“旋钮”。增大Q相对于R会使控制器更激进地跟踪参考但可能导致控制量饱和或震荡。增大R会使控制动作更平滑但响应变慢。通常从R的一个较小值开始逐渐增加Q直到获得满意的响应速度如果出现震荡再微调R。调整预测时域Hp和控制时域HcHp太短控制器“短视”可能为了短期利益如快速下降而采取损害长期稳定性的动作导致震荡。增加Hp通常能改善稳定性。Hc太短控制器“想象力”受限优化自由度小性能可能达不到最优。通常Hc设置为系统主要动态时间常数的2-3倍以采样周期计。引入输出权重或终端权重有时需要在预测时域末端施加一个终端代价P矩阵以保证稳定性。对于跟踪问题可以给输出误差加权。6.3 NMPC求解速度太慢无法实时运行问题fmincon一次优化耗时超过采样周期。优化策略热启动Warm Start如前所述这是必须的。用上一个周期的解作为当前周期的初始猜测。缩短时域在保证性能的前提下尽可能减少Hp和Hc。简化模型考虑使用降阶模型或拟线性化模型进行预测。使用更高效的求解器探索专用的实时NMPC求解器如acados、CasADi等它们比通用的fmincon快得多。显式MPC仅限线性对于线性系统带约束的MPC可以离线计算出状态分区和对应的控制律在线时只需查表速度极快。Matlab有explicitMPC功能。6.4 模型失配与鲁棒性问题实际系统与预测模型有差异导致控制性能下降甚至不稳定。增强鲁棒性增加输出反馈校正在预测时不是单纯用模型开环滚动而是用当前时刻的实际测量输出与模型预测输出的误差来修正未来的预测。这是MPC标准做法的一部分称为“反馈校正”或“摄动模型”。在手动实现时可以在每个周期初用测量值更新状态估计x(k)。使用鲁棒MPCRMPC在设计时考虑模型的不确定性范围求解一个 min-max 优化问题。这更复杂计算量更大。在线模型更新自适应MPC如果系统参数缓慢变化可以结合系统辨识方法在线更新预测模型。6.5 状态不可测与状态估计问题不是所有状态变量都能直接测量例如速度、某些内部温度。解决方案需要设计一个状态观测器。最常用的是卡尔曼滤波器对于线性系统或扩展卡尔曼滤波器EKF对于非线性系统。观测器根据可测量的输出y和控制输入u实时估计出全部状态x_hat然后将x_hat提供给MPC作为初始状态进行预测。在Matlab中对于线性系统可以使用kalman函数设计观测器。在Simulink中MPC Controller模块内部就集成了状态估计器默认使用稳态卡尔曼滤波器你只需要提供模型和测量噪声协方差。最后一个非常实用的建议是构建一个完整的、模块化的仿真测试平台。将你的MPC控制器、被控对象模型可以比预测模型更精细以测试鲁棒性、参考信号生成器、扰动和噪声模块都封装好。在投入实际应用前在这个平台上进行充分的测试包括标称性能测试、鲁棒性测试参数摄动、抗干扰测试等。这能帮你提前发现并解决大部分问题。MPC是一个强大的工具但它的有效性建立在准确的模型、合理的参数和严谨的工程实现之上。从简单的线性例子开始逐步增加复杂度是掌握它的最佳路径。