MPC模型预测控制：原理、设计与MATLAB实现

模型预测控制是一种先进的控制策略，特别适用于处理多变量约束控制问题。

1. MPC基本原理与核心思想

1.1 基本概念

MPC采用滚动优化和反馈校正的策略：

• 预测模型：使用系统数学模型预测未来动态
• 优化求解：在每个采样时刻求解有限时域优化问题
• 滚动实施：只实施当前时刻的控制量，下一时刻重新优化

1.2 MPC优势特点

• 直接处理多变量系统
• 显式处理约束（输入、输出、状态约束）
• 适用于非最小相位系统、大时滞系统
• 直观易懂，工程易接受

2. MPC数学 formulation

2.1 基本优化问题

\(min J = ∑‖y(k+i|k) - r(k+i)‖²_Q + ∑‖Δu(k+i|k)‖²_R\)

满足：

\(x(k+1) = Ax(k) + Bu(k)\)
\(y(k) = Cx(k)\)
\(u_min ≤ u(k) ≤ u_max\)
\(Δu_min ≤ Δu(k) ≤ Δu_max\)
\(y_min ≤ y(k) ≤ y_max\)

2.2 预测方程推导

对于状态空间模型：

% 预测方程构建
function [Y, U] = buildPredictionMatrices(A, B, C, Np, Nc)
    % Np: 预测时域, Nc: 控制时域
    
    n = size(A,1);  % 状态维度
    m = size(B,2);  % 输入维度
    p = size(C,1);  % 输出维度
    
    % 构建预测矩阵
    F = zeros(p*Np, n);
    Phi = zeros(p*Np, m*Nc);
    
    for i = 1:Np
        % 输出预测矩阵
        F((i-1)*p+1:i*p, :) = C * A^i;
        
        for j = 1:min(i,Nc)
            % 控制作用矩阵
            Phi((i-1)*p+1:i*p, (j-1)*m+1:j*m) = C * A^(i-j) * B;
        end
    end
end

3. MATLAB MPC工具箱使用

3.1 系统建模与MPC控制器创建

%% MPC控制器设计示例
clear; close all; clc;

% 1. 定义系统模型（离散状态空间）
Ts = 0.1; % 采样时间
A = [1, Ts; 0, 1];
B = [0.5*Ts^2; Ts];
C = [1, 0];
D = 0;
sys = ss(A, B, C, D, Ts);

% 2. 创建MPC控制器
mpcobj = mpc(sys, Ts);

% 3. 设置MPC参数
mpcobj.PredictionHorizon = 20;    % 预测时域
mpcobj.ControlHorizon = 3;        % 控制时域
mpcobj.Weights.OutputVariables = 1;   % 输出权重
mpcobj.Weights.ManipulatedVariablesRate = 0.1; % 控制变化权重

% 4. 设置约束
mpcobj.MV.Min = -2;      % 控制量下限
mpcobj.MV.Max = 2;       % 控制量上限
mpcobj.MV.RateMin = -1;  % 控制变化率下限
mpcobj.MV.RateMax = 1;   % 控制变化率上限
mpcobj.OV.Min = -10;     % 输出下限
mpcobj.OV.Max = 10;      % 输出上限

% 5. 显示控制器信息
display(mpcobj);

3.2 仿真验证

%% MPC闭环仿真
Tf = 10;        % 仿真时间
Tsteps = Tf/Ts; % 仿真步数

% 参考信号
r = 5 * ones(Tsteps, 1);
% 在中间时刻改变参考值
r(floor(Tsteps/2):end) = 2;

% 初始化
y = zeros(Tsteps, 1);
u = zeros(Tsteps, 1);
x = zeros(2, Tsteps);

% 创建仿真选项
options = mpcsimopt();
options.unmeas = [];     % 无不可测扰动
options.inputnoise = []; % 无输入噪声
options.outputnoise = 0.01 * randn(Tsteps,1); % 输出噪声

% 运行仿真
[Y, T, U, X] = sim(mpcobj, Tsteps, r, [], options);

% 绘制结果
figure('Position', [100, 100, 1200, 800]);

subplot(2,2,1);
plot(T, Y, 'b-', T, r, 'r--', 'LineWidth', 2);
xlabel('时间 (s)'); ylabel('输出');
title('系统输出 vs 参考信号');
legend('输出', '参考', 'Location', 'best');
grid on;

subplot(2,2,2);
plot(T, U, 'g-', 'LineWidth', 2);
xlabel('时间 (s)'); ylabel('控制量');
title('控制输入');
grid on;

subplot(2,2,3);
plot(T, X(:,1), 'm-', T, X(:,2), 'c-', 'LineWidth', 2);
xlabel('时间 (s)'); ylabel('状态');
title('系统状态');
legend('状态1', '状态2', 'Location', 'best');
grid on;

subplot(2,2,4);
error = r - Y;
plot(T, error, 'r-', 'LineWidth', 2);
xlabel('时间 (s)'); ylabel('误差');
title('跟踪误差');
grid on;

4. 手动实现MPC算法

4.1 QP问题 formulation

function [u_opt, cost] = mpcQP(A, B, C, x0, Np, Nc, Q, R, S, umin, umax, dumin, dumax, ymin, ymax, r)
    % MPC QP求解器
    % A,B,C: 系统矩阵
    % x0: 当前状态
    % Np: 预测时域, Nc: 控制时域
    % Q,R,S: 权重矩阵
    % 约束: umin, umax, dumin, dumax, ymin, ymax
    % r: 参考轨迹
    
    n = size(A,1); m = size(B,2); p = size(C,1);
    
    % 构建预测矩阵
    [F, Phi] = buildPredictionMatrices(A, B, C, Np, Nc);
    
    % QP问题: min 0.5*U'*H*U + f'*U
    H = Phi' * kron(eye(Np), Q) * Phi + kron(eye(Nc), R);
    f = (F*x0 - r)' * kron(eye(Np), Q) * Phi;
    
    % 控制增量约束: ΔU = U - U_prev
    A_delta = [];
    b_delta = [];
    for k = 1:Nc
        if k == 1
            A_delta = [A_delta; eye(m), zeros(m, (Nc-1)*m)];
            b_delta = [b_delta; dumax; -dumin];
        else
            temp = [zeros(m, (k-1)*m), eye(m), -eye(m), zeros(m, (Nc-k)*m)];
            A_delta = [A_delta; temp];
            b_delta = [b_delta; dumax; -dumin];
        end
    end
    
    % 输出约束
    A_y = Phi;
    b_y_max = ymax - F*x0;
    b_y_min = -(ymin - F*x0);
    
    % 控制量约束
    A_u = kron(eye(Nc), eye(m));
    b_u_max = umax;
    b_u_min = -umin;
    
    % 合并约束
    A_ineq = [A_delta; A_y; -A_y; A_u; -A_u];
    b_ineq = [b_delta; b_y_max; b_y_min; b_u_max; b_u_min];
    
    % 求解QP
    options = optimoptions('quadprog', 'Display', 'off');
    U_opt = quadprog(H, f', A_ineq, b_ineq, [], [], [], [], [], options);
    
    u_opt = U_opt(1:m);  % 只取第一个控制量
    cost = 0.5*U_opt'*H*U_opt + f*U_opt;
end

4.2 完整MPC仿真

%% 手动MPC实现 - 完整仿真
function manualMPCSimulation()
    % 系统参数
    Ts = 0.1;
    A = [0.8, 0.2; -0.1, 0.9];
    B = [0.5; 0.8];
    C = [1, 0];
    sys = ss(A, B, C, 0, Ts);
    
    % MPC参数
    Np = 15;    % 预测时域
    Nc = 3;     % 控制时域
    Q = 10;     % 输出权重
    R = 0.1;    % 控制权重
    
    % 约束
    umin = -2; umax = 2;
    dumin = -1; dumax = 1;
    ymin = -5; ymax = 5;
    
    % 仿真参数
    Tf = 10;
    Nsim = Tf/Ts;
    t = (0:Nsim-1)' * Ts;
    
    % 参考信号
    r = 3 * ones(Nsim, 1);
    r(floor(Nsim/3):floor(2*Nsim/3)) = 1;
    
    % 初始化
    x = zeros(2, Nsim);
    y = zeros(Nsim, 1);
    u = zeros(Nsim, 1);
    u_prev = 0;
    
    % MPC主循环
    for k = 1:Nsim-1
        % 获取当前状态（假设全状态可测）
        x_current = x(:, k);
        
        % 构建参考轨迹
        r_traj = r(k:min(k+Np-1, Nsim));
        if length(r_traj) < Np
            r_traj = [r_traj; r_traj(end)*ones(Np-length(r_traj),1)];
        end
        
        % 求解MPC问题
        [u_opt, cost] = mpcQP(A, B, C, x_current, Np, Nc, Q, R, 0, ...
                             umin, umax, dumin, dumax, ymin, ymax, r_traj);
        
        % 实施控制
        u(k) = u_opt;
        
        % 系统仿真
        x(:, k+1) = A * x(:, k) + B * u(k);
        y(k+1) = C * x(:, k+1);
        
        % 添加过程噪声（可选）
        % x(:, k+1) = x(:, k+1) + 0.01 * randn(2,1);
    end
    
    % 绘制结果
    plotMPCResults(t, y, u, r, x');
end

function plotMPCResults(t, y, u, r, x)
    figure('Position', [100, 100, 1200, 600]);
    
    subplot(2,2,1);
    plot(t, y, 'b-', t, r, 'r--', 'LineWidth', 2);
    xlabel('时间 (s)'); ylabel('输出');
    title('MPC控制: 输出跟踪');
    legend('输出', '参考', 'Location', 'best');
    grid on;
    
    subplot(2,2,2);
    plot(t, u, 'g-', 'LineWidth', 2);
    xlabel('时间 (s)'); ylabel('控制量');
    title('控制输入');
    grid on;
    
    subplot(2,2,3);
    plot(t, x(:,1), 'm-', t, x(:,2), 'c-', 'LineWidth', 2);
    xlabel('时间 (s)'); ylabel('状态');
    title('系统状态');
    legend('状态1', '状态2', 'Location', 'best');
    grid on;
    
    subplot(2,2,4);
    error = r - y;
    plot(t, error, 'r-', 'LineWidth', 1);
    xlabel('时间 (s)'); ylabel('误差');
    title('跟踪误差');
    grid on;
end

5. 非线性MPC实现

5.1 非线性系统MPC

%% 非线性MPC示例 - 倒立摆控制
function nonlinearMPCExample()
    % 倒立摆参数
    M = 1.0;    % 小车质量
    m = 0.1;    % 摆杆质量  
    l = 0.5;    % 摆杆长度
    g = 9.81;   % 重力加速度
    
    % 离散化参数
    Ts = 0.05;
    Tf = 5;
    Nsim = Tf/Ts;
    
    % MPC参数
    Np = 20;
    Nc = 3;
    Q = diag([10, 1, 10, 1]);  % 状态权重
    R = 0.1;                   % 控制权重
    
    % 初始状态 [x, x_dot, theta, theta_dot]
    x0 = [0; 0; pi-0.1; 0];  % 接近直立位置
    x = zeros(4, Nsim);
    x(:,1) = x0;
    
    % 参考状态（直立位置）
    x_ref = [0; 0; pi; 0];
    
    % 控制约束
    u_min = -20; u_max = 20;
    
    % 非线性MPC主循环
    options = optimoptions('fmincon', 'Display', 'off', ...
                          'Algorithm', 'sqp');
    
    for k = 1:Nsim-1
        % 当前状态
        x_current = x(:, k);
        
        % 求解非线性MPC
        u_opt = solveNonlinearMPC(@pendulumDynamics, x_current, x_ref, ...
                                 Np, Nc, Ts, Q, R, u_min, u_max, options);
        
        % 系统仿真（使用Runge-Kutta）
        [~, x_traj] = ode45(@(t,x) pendulumDynamics(x, u_opt), ...
                           [0, Ts], x_current);
        x(:, k+1) = x_traj(end, :)';
    end
    
    % 绘制结果
    plotPendulumResults((0:Nsim-1)*Ts, x, u_opt*ones(1,Nsim));
end

function dx = pendulumDynamics(x, u)
    % 倒立摆非线性动力学
    M = 1.0; m = 0.1; l = 0.5; g = 9.81;
    
    theta = x(3);
    dtheta = x(4);
    F = u;
    
    % 动力学方程
    denom = M + m*sin(theta)^2;
    
    dx1 = x(2);
    dx2 = (F + m*sin(theta)*(l*dtheta^2 + g*cos(theta))) / denom;
    dx3 = x(4);
    dx4 = (-F*cos(theta) - m*l*dtheta^2*cos(theta)*sin(theta) - (M+m)*g*sin(theta)) / (l*denom);
    
    dx = [dx1; dx2; dx3; dx4];
end

function u_opt = solveNonlinearMPC(dynamics, x0, x_ref, Np, Nc, Ts, Q, R, u_min, u_max, options)
    % 非线性MPC求解器
    
    % 初始猜测
    U0 = zeros(Nc, 1);
    
    % 优化问题
    costFunc = @(U) nonlinearCost(U, dynamics, x0, x_ref, Np, Nc, Ts, Q, R);
    
    % 约束
    A = []; b = []; Aeq = []; beq = [];
    lb = u_min * ones(Nc, 1);
    ub = u_max * ones(Nc, 1);
    
    % 求解
    U_opt = fmincon(costFunc, U0, A, b, Aeq, beq, lb, ub, [], options);
    u_opt = U_opt(1);
end

function J = nonlinearCost(U, dynamics, x0, x_ref, Np, Nc, Ts, Q, R)
    % 非线性MPC代价函数
    J = 0;
    x = x0;
    
    for k = 1:Np
        if k <= Nc
            u = U(k);
        else
            u = U(end);  % 控制时域外保持最后值
        end
        
        % 非线性仿真一步
        [~, x_traj] = ode45(@(t,x) dynamics(x, u), [0, Ts], x);
        x = x_traj(end, :)';
        
        % 累积代价
        J = J + (x - x_ref)' * Q * (x - x_ref) + u' * R * u;
    end
end

6. 实际应用考虑

6.1 抗扰与鲁棒性

%% 鲁棒MPC设计考虑
function robustMPCDesign()
    % 1. 干扰估计与补偿
    % 使用扰动观测器
    function d_est = disturbanceObserver(x, u, x_prev, A, B, L)
        x_pred = A * x_prev + B * u;
        d_est = L * (x - x_pred);
    end
    
    % 2. 软约束处理
    function cost = softConstraintsCost(epsilon, weight)
        % 对输出约束添加松弛变量
        cost = weight * sum(epsilon.^2);
    end
    
    % 3. 终端约束与代价
    function addTerminalConstraints()
        % 保证闭环稳定性
        % 方法: 终端代价函数、终端约束集
    end
end

6.2 实时实现考虑

%% 实时MPC实现技巧
function realTimeMPCTips()
    % 1. 显式MPC（离线计算）
    % 将MPC控制律预先计算为分段仿射函数
    
    % 2. 热启动优化
    % 使用上一时刻解作为初始猜测
    
    % 3. 简化模型
    % 使用降阶模型或线性参数变化模型
    
    % 4. 异步MPC
    % 根据系统动态调整MPC更新频率
end

7. MPC性能分析

7.1 稳定性与鲁棒性分析

%% MPC性能分析工具
function MPCAnalysis()
    % 1. 闭环极点分析
    % 计算MPC控制下的系统极点
    
    % 2. 灵敏度分析
    % 分析模型失配对性能的影响
    
    % 3. 约束激活分析
    % 统计约束激活频率和程度
    
    % 4. 计算负荷分析
    % 评估实时计算需求
end

参考代码 MPC模型预测控制 www.youwenfan.com/contentcnj/65398.html

总结

MPC的核心优势在于其处理多变量约束控制问题的能力。通过MATLAB实现时：

关键设计步骤：

1. 系统建模：获得准确的预测模型
2. 参数整定：调整预测时域、控制时域和权重
3. 约束设置：合理设置物理约束
4. 实时优化：选择高效的优化算法

MATLAB工具选择：

• MPC工具箱：适合快速原型开发和工业应用
• 手动实现：适合研究、教育和特殊需求
• 非线性MPC：适合复杂非线性系统

实际应用提示：

• 从简单模型开始，逐步增加复杂性
• 仔细验证约束的可行性
• 考虑计算实时性要求
• 加入适当的鲁棒性措施