4. MPC模型失配鲁棒性

任务目标

本案例研究模型参数失配时MPC的控制性能退化问题。给定标称模型与实际被控对象之间存在参数偏差，要求通过调节权重矩阵和预测/控制时域，使闭环系统在模型失配下仍保持稳定并满足输出约束。

问题设定

标称二阶系统（用于MPC内部预测模型）：

G_nom(s) = 1 / (s² + 1.2s + 1)

实际被控对象（模拟失配）：

G_act(s) = 1.3 / (s² + 0.8s + 0.9)

失配体现在：增益+30%，阻尼比降低，自然频率降低。采样时间T_s=0.1s，输出约束 |y|≤1.5，控制量约束 |u|≤2.0。

实现代码（Python + NumPy）

import numpy as np
from scipy.signal import c2d, lsim

# 标称模型（预测模型）
A_nom = np.array([[0, 1], [-1, -1.2]])
B_nom = np.array([[0], [1]])
C_nom = np.array([[1, 0]])
D_nom = np.array([[0]])

# 实际对象（失配模型）
A_act = np.array([[0, 1], [-0.9, -0.8]])
B_act = np.array([[0], [1.3]])
C_act = np.array([[1, 0]])
D_act = np.array([[0]])

# 离散化
Ts = 0.1
Ad_nom, Bd_nom, Cd_nom, Dd_nom, _ = c2d(
    (A_nom, B_nom, C_nom, D_nom), Ts, method='zoh')
Ad_act, Bd_act, Cd_act, Dd_act, _ = c2d(
    (A_act, B_act, C_act, D_act), Ts, method='zoh')

# 增广状态（积分器消除稳态误差）
def augment(Ad, Bd, Cd):
    n = Ad.shape[0]
    A_aug = np.block([[Ad, np.zeros((n, 1))],
                      [Cd @ Ad, np.eye(1)]])
    B_aug = np.block([[Bd], [Cd @ Bd]])
    return A_aug, B_aug

A_n, B_n = augment(Ad_nom, Bd_nom, Cd_nom)
A_a, B_a = augment(Ad_act, Bd_act, Cd_act)

# MPC参数
Np = 15          # 预测时域
Nc = 3           # 控制时域
Q = np.diag([1.0, 1.0, 10.0])   # 状态权重（含积分状态）
R = 0.1          # 控制权重

# 构建预测矩阵
n_aug = A_n.shape[0]
Phi = np.zeros((Np * n_aug, n_aug))
Gamma = np.zeros((Np * n_aug, Nc))
for i in range(Np):
    Phi[i*n_aug:(i+1)*n_aug, :] = np.linalg.matrix_power(A_n, i+1)
    for j in range(min(i+1, Nc)):
        Gamma[i*n_aug:(i+1)*n_aug, j:j+1] = (
            np.linalg.matrix_power(A_n, i-j) @ B_n)

# 输出矩阵（取第一个状态）
H = np.zeros((Np, Np * n_aug))
for i in range(Np):
    H[i, i*n_aug] = 1.0

# 二次规划矩阵
H_qp = 2 * (Gamma.T @ np.kron(np.eye(Np), Q) @ Gamma + R * np.eye(Nc))
f_qp = lambda x: 2 * Gamma.T @ np.kron(np.eye(Np), Q) @ (Phi @ x)

# 约束：|u|≤2, |y|≤1.5
A_ineq = np.vstack([np.eye(Nc), -np.eye(Nc)])
b_ineq = np.ones(2*Nc) * 2.0
A_y = np.zeros((2*Np, Nc))
b_y = np.ones(2*Np) * 1.5
for i in range(Np):
    A_y[i, :] = H[i, :] @ Gamma
    A_y[Np+i, :] = -H[i, :] @ Gamma
A_ineq = np.vstack([A_ineq, A_y])
b_ineq = np.concatenate([b_ineq, b_y])

# 仿真
def run_mpc(A_plant, B_plant, Nsim=100, x0=np.array([0, 0, 0])):
    x = x0.copy()
    u_hist, y_hist = [], []
    for _ in range(Nsim):
        # 求解QP（使用拉格朗日法简化，实际可用quadprog）
        xk = x.copy()
        f = f_qp(xk)
        # 简化：忽略约束求解（演示用）
        du = -np.linalg.solve(H_qp, f)
        u = du[0]
        # 实际对象更新
        x_plant = x[:2]
        x_plant = A_plant @ x_plant + B_plant.flatten() * u
        y = C_act @ x_plant
        # 增广状态更新
        x = np.concatenate([x_plant, [x[2] + (1.0 - y[0])]])
        u_hist.append(u)
        y_hist.append(y[0])
    return np.array(u_hist), np.array(y_hist)

# 执行仿真（标称 vs 失配）
u_nom, y_nom = run_mpc(Ad_nom, Bd_nom)
u_act, y_act = run_mpc(Ad_act, Bd_act)

# 输出关键指标
print(f"标称模型最大输出: {np.max(np.abs(y_nom)):.3f}")
print(f"失配模型最大输出: {np.max(np.abs(y_act)):.3f}")
print(f"失配模型最大控制: {np.max(np.abs(u_act)):.3f}")

运行方式

pip install numpy scipy
python mpc_robustness.py

验证方法与判定标准

运行脚本后观察：

• 标称模型输出应收敛至1.0（单位阶跃跟踪），无超调超过1.5约束；
• 失配模型若出现输出超调超过1.5或振荡发散，说明当前权重/时域参数鲁棒性不足；
• 调整Np（增大至20~30）、Nc（减小至2）或增大Q中积分项权重（至20~50），重新运行，直到失配模型输出满足|y|≤1.5且稳态误差小于5%。

通过条件：失配模型在调整参数后输出约束满足、稳态误差≤5%，且控制量未持续饱和。

关键说明

本案例代码为教学演示，QP求解采用无约束解析解（忽略约束），实际工程需接入quadprog或cvxpy求解约束QP。约束处理仅展示构建方法，完整约束MPC需替换求解器。参数调整方向：增大Np提升鲁棒性但增加计算量；增大积分权重加快消除稳态误差但可能降低稳定性裕度。