滚动时域控制（MPC）MATLAB 代码实例

线性 MPC 示例，对象为 双积分器系统（位置 + 速度），带输入约束。代码包含：

• 离散状态空间模型
• 预测模型构造（提升矩阵）
• 二次规划求解（quadprog）
• 滚动时域仿真循环
• 结果可视化

完整代码

%% 滚动时域控制（MPC）示例 —— 双积分器系统
clc; clear; close all;

%% 1. 系统模型（连续 → 离散）
% 双积分器：dx/dt = [0 1; 0 0] x + [0; 1] u
% 采样时间 Ts
Ts = 0.1;
Ac = [0 1; 0 0];
Bc = [0; 1];
Cc = eye(2);  % 输出为全状态
% 离散化（零阶保持）
sys_c = ss(Ac, Bc, Cc, []);
sys_d = c2d(sys_c, Ts, 'zoh');
[A, B, ~, ~] = ssdata(sys_d);

%% 2. MPC 参数
Np = 10;            % 预测时域
Nc = 5;             % 控制时域（通常 ≤ Np）
nx = size(A,1);     % 状态维数 = 2
nu = size(B,2);     % 输入维数 = 1

% 权重矩阵
Q = diag([10, 1]);   % 状态代价（位置权重高）
R = 0.1;             % 输入代价

% 约束
umin = -1;  umax = 1;   % 输入幅值约束
xmin = [-inf; -inf]; xmax = [inf; inf];  % 状态无约束（可加）

%% 3. 构建预测矩阵（提升形式）
% 预测模型：X_pred = F * x0 + G * U
% 其中 X_pred = [x(k+1); x(k+2); ...; x(k+Np)]
%       U = [u(k); u(k+1); ...; u(k+Nc-1)]

% F 矩阵：Np*nx × nx
F = zeros(Np*nx, nx);
% G 矩阵：Np*nx × Nc*nu
G = zeros(Np*nx, Nc*nu);

% 临时变量
temp_F = eye(nx);
temp_G = zeros(nx, Nc*nu);
for i = 1:Np
    % 第 i 步预测
    temp_F = A * temp_F;   % A^i
    F((i-1)*nx+1:i*nx, :) = temp_F;
    
    % 构造 G 的第 i 块行
    % 第 j 个控制量对第 i 步的影响：A^(i-j)*B  (j=1..min(i,Nc))
    for j = 1:min(i, Nc)
        G((i-1)*nx+1:i*nx, (j-1)*nu+1:j*nu) = ...
            A^(i-j) * B;
    end
end

%% 4. 构造二次规划矩阵
% 目标函数：J = (X_ref - X_pred)' * Q_bar * (X_ref - X_pred) + U' * R_bar * U
% 其中 Q_bar = blkdiag(Q, Q, ..., Q) (Np个)
%       R_bar = blkdiag(R, R, ..., R) (Nc个)
Q_bar = kron(eye(Np), Q);
R_bar = kron(eye(Nc), R);

% 将目标函数化为标准二次型：J = 0.5 * U' * H * U + f' * U + const
% 代入 X_pred = F*x0 + G*U，忽略常数项：
% J = (F*x0 + G*U)' * Q_bar * (F*x0 + G*U) + U' * R_bar * U
%   = U' * (G'*Q_bar*G + R_bar) * U + 2 * (F*x0)' * Q_bar * G * U + const
H = 2 * (G' * Q_bar * G + R_bar);   % 注意 quadprog 使用 0.5*U'*H*U，所以这里乘2
% 实际上 quadprog 要求 H 为二次项系数的一半，所以 H = 2*(G'*Q_bar*G + R_bar) 才是正确的
% 更严谨：quadprog 的 H 应满足 0.5*x'*H*x，所以 H = 2*(G'*Q_bar*G + R_bar)
% 但为清晰，我们直接使用 quadprog 的标准形式：min 0.5*x'*H*x + f'*x
% 故 H = 2*(G'*Q_bar*G + R_bar)

% 参考状态（设为原点）
xref = zeros(Np*nx, 1);
% 线性项 f = -2 * (xref - F*x0)' * Q_bar * G = 2 * (F*x0)' * Q_bar * G (因为 xref=0)
% 写成 f'*U 的形式，f = -2 * G' * Q_bar * (xref - F*x0) = 2 * G' * Q_bar * F * x0

%% 5. 仿真参数
Tsim = 50;          % 仿真步数
x0 = [2; 0];        % 初始状态 [位置; 速度]
x = x0;             % 当前状态
u_hist = zeros(Tsim, 1);
x_hist = zeros(Tsim, nx);
x_hist(1,:) = x';

%% 6. MPC 主循环（滚动时域）
options = optimoptions('quadprog', 'Display', 'off', 'Algorithm', 'interior-point-convex');

for k = 1:Tsim
    % 当前状态
    xk = x;
    
    % 构造线性项 f (依赖于当前状态)
    f = 2 * G' * Q_bar * F * xk;
    
    % 输入约束：umin <= u <= umax
    lb = umin * ones(Nc, 1);
    ub = umax * ones(Nc, 1);
    
    % 状态约束（如有需要，可在此添加）
    % 无状态约束时省略 Aineq, bineq
    
    % 求解二次规划
    U_opt = quadprog(H, f, [], [], [], [], lb, ub, [], options);
    
    % 取第一个控制量作用于系统
    u = U_opt(1);
    u_hist(k) = u;
    
    % 更新系统状态（真实系统）
    x = A * x + B * u;
    if k < Tsim
        x_hist(k+1,:) = x';
    end
end

%% 7. 结果可视化
figure;
subplot(2,1,1);
stairs(0:Tsim-1, u_hist, 'LineWidth', 1.5); hold on;
yline([umin, umax], 'r--', 'LineWidth', 1);
xlabel('时间步 k'); ylabel('控制输入 u');
title('MPC 控制输入'); legend('u', '约束'); grid on;

subplot(2,1,2);
plot(0:Tsim-1, x_hist(:,1), 'b-', 'LineWidth', 1.5); hold on;
plot(0:Tsim-1, x_hist(:,2), 'r--', 'LineWidth', 1.5);
xlabel('时间步 k'); ylabel('状态');
title('系统状态响应'); legend('位置', '速度'); grid on;

sgtitle('滚动时域控制 (MPC) 仿真 —— 双积分器');

代码说明

1. 系统模型

• 连续双积分器：$\ddot{x}=u$，状态 $x=[p,v{]}^T$
• 零阶保持离散化，采样时间 0.1s

2. MPC 参数

• 预测时域 $N_p=10$，控制时域 $N_c=5$
• 权重 $Q=\text{diag}(10,1)$（强调位置跟踪），$R=0.1$
• 输入约束 $u\in [-1,1]$

3. 预测矩阵

• 提升形式：$X_{\text{pred}}=F{x}_k+GU$
• 通过循环构造 $F$ 和 $G$ 矩阵

4. 二次规划

• 目标函数转化为标准 QP：${\min }_U\frac{1}{2}{U}^{T}HU+f^{T}U$
• 使用 quadprog 求解，只施加输入上下界约束

5. 滚动仿真

• 每个时刻求解 QP，取第一个控制量 $u(k)$ 作用于系统
• 更新状态，进入下一时刻

运行结果

• 控制输入：始终保持在 $[-1,1]$ 范围内
• 状态响应：位置从 2 逐渐回到 0，速度先增加后减小，系统稳定

参考代码 滚动时域控制的代码实例 www.youwenfan.com/contentcsv/81650.html

扩展建议