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

任务目标

实现一阶惯性系统（G(s) = K/(Ts+1)）的模型预测控制器，使其输出跟踪阶跃参考信号。通过仿真验证MPC的预测控制能力，并分析闭环响应特性。

系统模型与离散化

连续一阶惯性系统：

G(s) = K / (Ts + 1)

采用零阶保持器离散化（采样时间 Ts_sample），得到离散状态空间模型：

x(k+1) = a·x(k) + b·u(k)
y(k)   = c·x(k)

其中：a = exp(-Ts_sample/T)，b = K·(1-a)，c = 1。

MPC控制器设计

预测时域 Np，控制时域 Nc，目标函数：

J = Σ(y(k+i|k) - r(k+i))² + λ·Σ Δu(k+i)²

约束：控制增量 |Δu| ≤ Δu_max，控制量 |u| ≤ u_max。

完整实现代码

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

# 系统参数
K = 2.0          # 增益
T = 5.0          # 时间常数
Ts = 0.5         # 采样时间
a = np.exp(-Ts/T)
b = K * (1 - a)

# MPC参数
Np = 10          # 预测时域
Nc = 3           # 控制时域
lambda_u = 0.1   # 控制增量权重
du_max = 0.5     # 控制增量约束
u_max = 3.0      # 控制量约束

# 仿真参数
T_sim = 30.0
N_steps = int(T_sim / Ts)
r = 1.5          # 阶跃参考值

# 初始化
x = 0.0
u_prev = 0.0
y_history = []
u_history = []
r_history = []

def predict(x0, du_seq):
    """基于当前状态和控制增量序列预测未来输出"""
    x_pred = x0
    u_pred = u_prev
    y_pred = []
    for i in range(Np):
        if i < len(du_seq):
            u_pred += du_seq[i]
        u_pred = np.clip(u_pred, -u_max, u_max)
        x_pred = a * x_pred + b * u_pred
        y_pred.append(x_pred)
    return np.array(y_pred)

def cost(du_seq):
    """MPC目标函数"""
    y_pred = predict(x, du_seq)
    r_seq = np.ones(Np) * r
    du_seq_full = np.zeros(Np)
    du_seq_full[:Nc] = du_seq
    cost_track = np.sum((y_pred - r_seq)**2)
    cost_du = lambda_u * np.sum(du_seq_full**2)
    return cost_track + cost_du

# 仿真主循环
for k in range(N_steps):
    # 求解MPC优化问题
    du0 = np.zeros(Nc)
    bounds = [(-du_max, du_max)] * Nc
    res = minimize(cost, du0, method='SLSQP', bounds=bounds)
    du_opt = res.x
    
    # 应用第一个控制增量
    u_new = np.clip(u_prev + du_opt[0], -u_max, u_max)
    
    # 更新系统状态
    x = a * x + b * u_new
    
    # 记录数据
    y_history.append(x)
    u_history.append(u_new)
    r_history.append(r)
    u_prev = u_new

# 绘制结果
t = np.arange(N_steps) * Ts
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(t, y_history, 'b-', label='系统输出 y')
plt.plot(t, r_history, 'r--', label='参考信号 r')
plt.ylabel('输出 y')
plt.legend()
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(t, u_history, 'g-', label='控制量 u')
plt.ylabel('控制量 u')
plt.xlabel('时间 (s)')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

# 计算性能指标
y_arr = np.array(y_history)
steady_idx = int(N_steps * 0.7)
overshoot = (np.max(y_arr) - r) / r * 100
settling_idx = np.where(np.abs(y_arr - r) < 0.05 * r)[0]
settling_time = (settling_idx[-1] - settling_idx[0]) * Ts if len(settling_idx) > 0 else None

print(f"超调量: {overshoot:.2f}%")
print(f"稳态误差: {abs(y_arr[-1] - r):.4f}")
if settling_time:
    print(f"调节时间(±5%): {settling_time:.1f}s")

运行方式

1. 安装依赖：pip install numpy scipy matplotlib
2. 将代码保存为 mpc_first_order.py 并运行：python mpc_first_order.py

验证方法与判定标准

• 跟踪性能：系统输出应平滑跟踪阶跃参考值 r=1.5，稳态误差小于 0.05
• 约束满足：控制量始终满足 |u| ≤ 3.0，控制增量满足 |Δu| ≤ 0.5
• 动态响应：输出无剧烈振荡，超调量通常小于 10%（可通过调整 λ_u 控制）
• 预期输出：程序运行后显示超调量约 0%（λ_u=0.1 时），稳态误差接近 0，调节时间约 8~12s

关键参数影响分析

常见问题

• 求解失败：若 SLSQP 求解器报错，可尝试增大 maxiter 参数或改用 trust-constr 方法
• 响应振荡：增大 λ_u 或减小 Nc 可抑制振荡
• 稳态误差不为零：检查系统模型是否准确，或考虑加入积分项（增量式MPC已包含积分作用）