企业官网建设流程全解析

用MATLAB手把手实现一个简易的MPC控制器（附完整代码与避坑指南）

在工业控制和自动化领域，模型预测控制（MPC）因其出色的多变量处理能力和约束处理特性而备受青睐。但对于初学者来说，从理论到实践的跨越往往充满挑战——复杂的矩阵运算、晦涩的优化问题，还有那些让人头疼的维度匹配错误。本文将用最直观的方式，带你用MATLAB从零搭建一个可运行的MPC控制器，避开那些教科书不会告诉你的"坑"。

1. 基础准备：理解MPC的核心骨架

MPC本质上是一个"预测-优化-执行"的循环过程。想象你在驾驶汽车时，会不断预测前方路况，调整方向盘和油门，这就是MPC的直观体现。数学上，这个循环包含三个关键环节：

• 预测模型：用状态空间方程描述系统动态
• 滚动优化：在每个时刻求解最优控制序列
• 反馈校正：用实际测量值修正预测误差

对于离散线性系统，状态空间方程可以表示为：

% 离散状态空间方程示例 Ad = [0.8 0.2; -0.1 0.9]; % 状态矩阵 Bd = [0.5; 0.3]; % 输入矩阵 Cd = [1 0]; % 输出矩阵

提示：实际工程中，通常先用c2d函数将连续模型离散化，采样时间的选择会影响控制性能

2. 预测矩阵构建：从数学到代码

预测时域内的系统行为可以表示为：

X = F * x0 + Φ * U + E

其中F和Φ是最关键的预测矩阵。让我们用MATLAB构建它们：

function [F, Phi, E] = build_prediction_matrices(Ad, Bd, Cd, Dd, Np, Nc) nx = size(Ad,1); % 状态维度 nu = size(Bd,2); % 输入维度 ny = size(Cd,1); % 输出维度 F = zeros(Np*ny, nx); Phi = zeros(Np*ny, Nc*nu); E = zeros(Np*ny, size(Dd,2)); % 构建F矩阵 F(1:ny,:) = Cd * Ad; for k = 2:Np F((k-1)*ny+1:k*ny,:) = F((k-2)*ny+1:(k-1)*ny,:) * Ad; end % 构建Φ矩阵 for i = 1:Np for j = 1:min(i,Nc) Phi((i-1)*ny+1:i*ny, (j-1)*nu+1:j*nu) = ... Cd * Ad^(i-j) * Bd; end end % 构建扰动矩阵E（假设有常数扰动） E(1:ny,:) = Cd * Dd; for k = 2:Np E((k-1)*ny+1:k*ny,:) = Cd * Ad^(k-1) * Dd + E((k-2)*ny+1:(k-1)*ny,:); end end

常见错误：矩阵维度不匹配通常源于Ad,Bd,Cd维度定义不一致，建议在代码开头添加维度检查

3. 优化问题求解：QP问题的实战转换

MPC的核心是将控制问题转化为二次规划(QP)：

min J = U'*H*U + 2*U'*G s.t. U_min ≤ U ≤ U_max

对应的MATLAB实现：

function [U_opt, flag] = solve_mpc_qp(H, G, U_min, U_max) nu = length(U_min); Nc = nu / size(U_min,1); % 构建不等式约束 A*U ≤ b A = [eye(nu); -eye(nu)]; b = [U_max; -U_min]; % 使用quadprog求解 options = optimoptions('quadprog', 'Display', 'none'); [U_opt, ~, exitflag] = quadprog(H, G, A, b, [], [], [], [], [], options); flag = exitflag == 1; if ~flag warning('QP求解失败，尝试无约束解'); U_opt = -H\G; % 退化为无约束解 end end

参数选择经验值：

4. 闭环实现与调试技巧

完整的MPC闭环实现需要处理几个关键点：

function [u, x_pred] = mpc_controller(x, u_prev, Ad, Bd, Cd, Q, R, Np, Nc, umax, umin) % 构建预测矩阵 [F, Phi, E] = build_prediction_matrices(Ad, Bd, Cd, 0, Np, Nc); % 构建QP参数 Qe = kron(eye(Np), Q); Re = kron(eye(Nc), R); H = Phi'*Qe*Phi + Re; G = Phi'*Qe*(F*x); % 求解QP U_min = kron(ones(Nc,1), umin); U_max = kron(ones(Nc,1), umax); [U_opt, ~] = solve_mpc_qp(H, G, U_min, U_max); % 取第一个控制量 u = U_opt(1:length(umin)); x_pred = F*x + Phi*U_opt; end

调试中常见问题及解决方案：

1. QP求解失败：

   • 检查H矩阵的正定性
   • 尝试减小Nc或放松约束

2. 控制发散：

   • 增加Q权重强化状态调节
   • 检查模型(Ad,Bd)是否正确

3. 计算延迟：

   • 减少Np/Nc
   • 使用热启动技巧（用上一时刻的解作为初始值）

% 典型的主循环结构 for k = 1:sim_steps [u, x_pred] = mpc_controller(x, u, Ad, Bd, Cd, Q, R, Np, Nc, umax, umin); x = Ad*x + Bd*u; % 系统动态 y = Cd*x; % 记录数据 U_history(:,k) = u; X_history(:,k) = x; end

5. 进阶技巧：提升MPC性能的实用方法

当基础MPC运行稳定后，可以考虑以下优化：

1. 终端代价与约束：

% 在H矩阵中添加终端代价 P = dare(Ad, Bd, Q, R); % 求解Riccati方程 H = Phi'*Qe*Phi + Re + Phi(end-ny+1:end,:)'*P*Phi(end-ny+1:end,:);

2. 扰动估计与抗偏移：

% 增广状态空间包含扰动估计 Ad_aug = [Ad Bd; zeros(nu,nx) eye(nu)]; Bd_aug = [Bd; eye(nu)];

3. 软约束处理：

% 在QP中引入松弛变量 H = blkdiag(Phi'*Qe*Phi + Re, 1e6*eye(ny)); G = [Phi'*Qe*F*x; zeros(ny,1)]; A = [A zeros(size(A,1),ny); -Cd*Phi -eye(ny)]; b = [b; -ymin + Cd*F*x];

6. 完整代码实现与测试案例

以下是一个完整的SISO系统MPC控制示例：

%% 系统定义 sys = ss(tf([1],[1 1.5 1])); % 二阶振荡系统 Ts = 0.1; % 采样时间 sysd = c2d(sys, Ts); % 离散化 Ad = sysd.A; Bd = sysd.B; Cd = sysd.C; %% MPC参数 Np = 15; % 预测时域 Nc = 3; % 控制时域 Q = 10; % 输出权重 R = 0.1; % 输入权重 umax = 1; umin = -1; % 输入约束 %% 仿真设置 sim_time = 5; % 秒 sim_steps = sim_time/Ts; ref = 1; % 参考值 %% 初始化 x = zeros(size(Ad,1),1); u = 0; X_history = zeros(size(Ad,1), sim_steps); U_history = zeros(1, sim_steps); Y_history = zeros(1, sim_steps); %% 主循环 for k = 1:sim_steps % 调用MPC控制器 [u, ~] = mpc_controller(x, u, Ad, Bd, Cd, Q, R, Np, Nc, umax, umin); % 系统动态 x = Ad*x + Bd*u; y = Cd*x; % 记录数据 U_history(:,k) = u; X_history(:,k) = x; Y_history(:,k) = y; end %% 可视化 t = (0:sim_steps-1)*Ts; figure; subplot(2,1,1); plot(t,Y_history); ylabel('Output'); subplot(2,1,2); plot(t,U_history); ylabel('Control'); xlabel('Time(s)');

这个案例展示了如何控制一个二阶振荡系统跟踪阶跃信号。在实际项目中，你可能需要：

1. 添加参考跟踪项
2. 实现状态估计器（如Kalman滤波）
3. 处理输入速率约束
4. 加入积分动作消除稳态误差

7. 性能优化与实时实现

当需要在嵌入式设备上部署MPC时，考虑以下优化策略：

代码生成友好实现：

function [u, status] = mpc_autocode(x, u_prev, params) % 将矩阵计算展开为标量运算 % 适用于代码生成 persistent H cholH % 复用矩阵分解 if isempty(H) % 离线计算不变部分 [~, Phi, ~] = build_prediction_matrices(... params.Ad, params.Bd, params.Cd, 0, params.Np, params.Nc); H = Phi'*kron(eye(params.Np),params.Q)*Phi + kron(eye(params.Nc),params.R); cholH = chol(H,'lower'); end % 在线计算 Fx = compute_Fx(x, params); % 展开的矩阵乘法 G = Phi'*kron(eye(params.Np),params.Q)*Fx; % 求解三角系统 z = linsolve(cholH, G, struct('LT',true)); U_opt = -linsolve(cholH', z, struct('UT',true)); u = U_opt(1) + params.Kf*u_prev; % 可能添加滤波 status = 1; end

计算耗时优化技巧：

• 预计算不变矩阵（如H矩阵的Cholesky分解）
• 使用定点算术减少计算量
• 采用显式MPC（预先计算最优控制律的分区表达式）

在实时性要求高的场合，实测表明对于4状态2输入的系统，在STM32H7上单次QP求解可控制在1ms以内。