4. 非最小相位系统MPC控制

任务目标

针对具有右半平面零点(非最小相位)的连续系统,设计离散MPC控制器,实现输出跟踪,并分析零点位置对控制性能(超调、响应速度)的约束。

系统模型

选择典型非最小相位传递函数:

G(s) = (-s + 2) / (s² + s + 1)

该模型在s=2处有右半平面零点。采样时间Ts=0.1s,采用零阶保持器离散化。

MPC控制器设计

预测时域Np=10,控制时域Nc=3,输出权重Q=1,控制增量权重R=0.1。约束:|u|≤2,|Δu|≤0.5。

关键完整代码

% 非最小相位系统MPC控制
clear; clc;

% 系统定义
Ts = 0.1;
G = tf([-1 2], [1 1 1]);  % 非最小相位系统
Gd = c2d(G, Ts, 'zoh');
A = Gd.A; B = Gd.B; C = Gd.C; D = Gd.D;
[nx, nu] = size(B);

% MPC参数
Np = 10; Nc = 3;
Q = 1; R = 0.1;
u_min = -2; u_max = 2;
du_min = -0.5; du_max = 0.5;

% 构建预测矩阵
[Phi, Gamma] = build_prediction_matrices(A, B, C, Np, Nc);

% 构建QP矩阵
H = 2*(Gamma'*Q*Gamma + R*eye(Nc));
f = @(x, u_prev) 2*Gamma'*Q*(Phi*x - 0);  % 目标输出为0

% 约束矩阵
A_ineq = [tril(ones(Nc)); -tril(ones(Nc)); eye(Nc); -eye(Nc)];
b_ineq = @(u_prev) [du_max*ones(Nc,1) + u_prev*ones(Nc,1);
                    -du_min*ones(Nc,1) - u_prev*ones(Nc,1);
                    u_max*ones(Nc,1); -u_min*ones(Nc,1)];

% 仿真
T_sim = 10;
N_steps = T_sim / Ts;
x = zeros(nx, 1);
u = 0;
y_history = zeros(N_steps, 1);
u_history = zeros(N_steps, 1);

for k = 1:N_steps
    % 当前状态
    x_k = x;
    u_prev = u;
    
    % 求解QP
    H_eff = (H + H')/2;
    f_eff = f(x_k, u_prev);
    A_eff = A_ineq;
    b_eff = b_ineq(u_prev);
    
    options = optimoptions('quadprog', 'Display', 'off');
    [du_opt, ~] = quadprog(H_eff, f_eff, A_eff, b_eff, [], [], [], [], [], options);
    
    % 应用第一个控制增量
    du = du_opt(1);
    u = u_prev + du;
    u = max(u_min, min(u_max, u));
    
    % 系统更新
    x = A*x + B*u;
    y = C*x + D*u;
    
    % 记录
    y_history(k) = y;
    u_history(k) = u;
end

% 绘图
t = (0:N_steps-1)*Ts;
figure;
subplot(2,1,1);
plot(t, y_history, 'b-', 'LineWidth', 1.5);
ylabel('输出 y');
title('非最小相位系统MPC响应');
grid on;
subplot(2,1,2);
plot(t, u_history, 'r-', 'LineWidth', 1.5);
ylabel('控制 u');
xlabel('时间 (s)');
grid on;

% 辅助函数:构建预测矩阵
function [Phi, Gamma] = build_prediction_matrices(A, B, C, Np, Nc)
    [nx, ~] = size(A);
    Phi = zeros(Np, nx);
    Gamma = zeros(Np, Nc);
    
    % 构建Phi
    for i = 1:Np
        Phi(i, :) = C * A^i;
    end
    
    % 构建Gamma
    for i = 1:Np
        for j = 1:min(i, Nc)
            Gamma(i, j) = C * A^(i-j) * B;
        end
    end
end

运行方式

  1. 需要MATLAB环境(R2016b及以上),无需额外工具箱(quadprog来自Optimization Toolbox)
  2. 将代码保存为mpc_nonminimum.m,直接运行

验证方法与判定标准

预期输出特征:

  • 系统输出从初始值0开始,在阶跃参考(此处为0)激励下,输出呈现初始反向响应(先向参考反方向偏离,再折回),这是非最小相位系统的典型特征
  • 控制量u在约束范围内(|u|≤2),且变化率受限(|Δu|≤0.5)
  • 输出最终收敛到0(稳态无偏差)

性能限制分析:右半平面零点位于s=2,对应零点时间常数τz=0.5s。MPC预测时域Np=10×0.1=1s,大于2τz,控制器能够预见反向响应并提前调整控制量。若减小Np至小于0.5s,输出将出现明显超调甚至振荡,验证零点对响应速度的根本限制。

常见问题

  • quadprog报错:确认H矩阵对称正定,代码中已用(H+H')/2强制对称
  • 初始反向响应不明显:增大仿真时间或减小R权重,使控制器更积极
  • 约束冲突:若du_max过小导致无解,可适当放宽至0.8