📚 课程目录第1章一阶惯性系统MPC跟踪阶跃响应

1. 一阶惯性系统MPC跟踪阶跃响应

任务目标

为一阶惯性系统 G(s) = 1/(s+1) 设计模型预测控制器（MPC），使其输出跟踪阶跃参考信号，并对比不同预测时域 Np 下的跟踪性能差异。

系统离散化与MPC设计

采用零阶保持器以采样时间 Ts = 0.1s 对连续系统离散化，得到状态空间模型：

x(k+1) = Ad·x(k) + Bd·u(k)
y(k)   = Cd·x(k)

MPC在每个采样时刻求解有限时域优化问题：

min  Σ_{i=1}^{Np} (y(k+i|k) - r(k+i))² + λ·Σ_{i=0}^{Nc-1} Δu(k+i|k)²
s.t. u_min ≤ u ≤ u_max

其中 Np 为预测时域，Nc 为控制时域（取 Nc = Np），λ 为控制增量权重。通过滚动优化，仅执行当前时刻最优控制量 u(k)，下一时刻重新求解。

完整Python实现

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize

# 系统参数
Ts = 0.1
A = np.array([[-1.0]])
B = np.array([[1.0]])
C = np.array([[1.0]])
D = np.array([[0.0]])

# 离散化（零阶保持）
Ad = np.array([[np.exp(-Ts)]])
Bd = np.array([[1 - np.exp(-Ts)]])
Cd = C.copy()

# MPC参数
Np = 10          # 预测时域（可修改）
Nc = Np          # 控制时域
lam = 0.1        # 控制增量权重
u_min, u_max = -2.0, 2.0
du_max = 0.5     # 控制增量限幅

# 仿真参数
T_total = 5.0
N_steps = int(T_total / Ts)
r = 1.0          # 阶跃参考值
x0 = np.array([[0.0]])

def mpc_solve(xk, rk, Np, Nc):
    """求解MPC优化问题，返回当前控制增量du"""
    n = len(xk)
    # 构建预测矩阵
    A_pow = np.zeros((Np, n, n))
    A_pow[0] = Ad
    for i in range(1, Np):
        A_pow[i] = Ad @ A_pow[i-1]
    
    # 预测矩阵 F（自由响应）和 G（受控响应）
    F = np.zeros((Np, n))
    G = np.zeros((Np, Nc))
    for i in range(Np):
        F[i] = Cd @ A_pow[i]
        for j in range(min(i+1, Nc)):
            G[i, j] = Cd @ A_pow[i-j-1] @ Bd
    
    # 构建QP目标函数（转化为du的二次型）
    R = np.eye(Nc) * lam
    H = G.T @ G + R
    f = (F @ xk - rk).T @ G
    
    # 约束：du_min ≤ du ≤ du_max
    bounds = [(-du_max, du_max)] * Nc
    
    # 求解
    def cost(du):
        return 0.5 * du @ H @ du + f @ du
    
    res = minimize(cost, np.zeros(Nc), bounds=bounds, method='SLSQP')
    return res.x[0]

# 仿真主循环
x = x0.copy()
u_hist = []
y_hist = [float(Cd @ x)]
x_hist = [float(x[0])]

for k in range(N_steps):
    rk = np.array([[r]])
    du = mpc_solve(x, rk, Np, Nc)
    u_new = (u_hist[-1] if u_hist else 0) + du
    u_new = np.clip(u_new, u_min, u_max)
    
    # 更新状态
    x = Ad @ x + Bd * u_new
    u_hist.append(u_new)
    y_hist.append(float(Cd @ x))
    x_hist.append(float(x[0]))

# 绘制结果
t = np.arange(N_steps+1) * Ts
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(t, y_hist, 'b-', linewidth=2, label='MPC输出')
plt.plot(t, np.ones_like(t)*r, 'r--', label='参考值')
plt.xlabel('时间 (s)')
plt.ylabel('输出 y')
plt.title(f'Np={Np}')
plt.legend()
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(t[:-1], u_hist, 'g-', linewidth=2)
plt.xlabel('时间 (s)')
plt.ylabel('控制输入 u')
plt.title('控制输入')
plt.grid(True)
plt.tight_layout()
plt.show()

# 计算性能指标
y_arr = np.array(y_hist)
rise_time = None
settle_time = None
for i, y_val in enumerate(y_arr):
    if rise_time is None and y_val >= 0.9 * r:
        rise_time = i * Ts
    if settle_time is None and i > 5 and abs(y_val - r) < 0.02 * r:
        settle_time = i * Ts
        break

print(f"预测时域 Np = {Np}")
print(f"上升时间: {rise_time:.2f}s" if rise_time else "上升时间: 未达到")
print(f"调节时间(±2%): {settle_time:.2f}s" if settle_time else "调节时间: 未达到")
print(f"稳态误差: {abs(y_arr[-1] - r):.4f}")

运行与验证

运行方式：直接执行上述Python脚本。需要安装 numpy、scipy、matplotlib。

验证方法：

运行代码后观察输出曲线是否跟踪阶跃参考值1.0；
修改 Np 值（如分别设为5、10、20），重新运行并对比上升时间与调节时间；
预期输出：Np 较小时响应更快但可能振荡，Np 较大时响应更平缓但上升时间略长。稳态误差应趋近于0。

预测时域影响分析

将 Np 分别设为5、10、20运行脚本，可观察到：

Np=5：预测范围短，控制器更"短视"，响应速度快但超调量增大，控制输入变化更剧烈；
Np=10：平衡了快速性与稳定性，超调减小，调节时间适中；
Np=20：预测范围长，控制器更"保守"，响应平缓、无超调，但上升时间变长。

实际应用中需根据系统对快速性与平稳性的要求折中选择 Np。

常见问题

控制量饱和：若 u_min/u_max 设置过小，MPC输出可能长时间饱和，导致跟踪变慢。可增大限幅或调整 λ 权重。
数值问题：当 Np 很大时，H 矩阵可能接近奇异。可增加 λ 值改善条件数。
模型失配：实际系统参数与模型不一致时，可加入状态估计器或扰动补偿提升鲁棒性。