3. 滚转轴PID控制器响应速度优化

任务说明

本案例针对火箭滚转轴姿态控制回路，要求整定PID控制器参数，使系统对阶跃指令的响应满足：上升时间 t_r < 0.3 s，且超调量 M_p < 5%。被控对象为滚转通道简化模型，采用PD+积分结构，通过参数寻优获得满足双指标的组合。

被控对象与控制器结构

滚转轴动力学可近似为一阶惯性环节加执行机构延迟，传递函数为：

G(s) = 1 / (0.8s + 1) · e<sup>-0.05s</sup>

控制器采用并联PID形式：

u(t) = K<sub>p</sub>·e(t) + K<sub>i</sub>·∫e(τ)dτ + K<sub>d</sub>·de(t)/dt

其中 e(t) 为期望滚转角与当前滚转角之差。执行机构输出限幅为 ±30°。

参数整定实现

采用时域仿真结合网格搜索法，在 K_p∈[2,10]、K_i∈[0,5]、K_d∈[0.5,4] 范围内以步长0.1遍历，对每组参数执行阶跃仿真，筛选同时满足上升时间与超调量约束的参数组合，并从中选取超调量最小者作为最终参数。


import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# 被控对象参数
tau = 0.8      # 时间常数
delay = 0.05   # 纯延迟时间
sat_limit = 30.0  # 执行机构限幅

# 仿真参数
dt = 0.001
t_end = 3.0
t = np.arange(0, t_end, dt)
n_delay = int(delay / dt)

# 阶跃目标
setpoint = 10.0  # 期望滚转角（度）

def simulate_pid(kp, ki, kd):
    """对给定PID参数执行阶跃响应仿真，返回时间序列与输出"""
    n = len(t)
    y = np.zeros(n)          # 实际输出
    u = np.zeros(n)          # 控制器输出（限幅前）
    u_sat = np.zeros(n)      # 限幅后输出
    e_integral = 0.0
    e_prev = 0.0
    y_delay = np.zeros(n_delay + 1)  # 延迟缓冲

    for i in range(n):
        # 当前实际输出（考虑延迟）
        if i >= n_delay:
            y_actual = y_delay[i - n_delay]
        else:
            y_actual = 0.0

        # 误差
        e = setpoint - y_actual
        e_integral += e * dt
        # 微分项（对测量值微分，避免微分冲击）
        dedt = (y_actual - y_delay[i - 1]) / dt if i > 0 else 0.0

        # PID控制律
        u[i] = kp * e + ki * e_integral - kd * dedt

        # 执行机构限幅
        u_sat[i] = np.clip(u[i], -sat_limit, sat_limit)

        # 更新被控对象状态（一阶惯性）
        if i == 0:
            y[i] = 0.0
        else:
            # 使用前向欧拉法
            y[i] = y[i-1] + dt * (u_sat[i-1] - y[i-1]) / tau

        # 更新延迟缓冲
        y_delay[i] = y[i]

    return t, y

def evaluate_response(t, y):
    """计算上升时间与超调量"""
    # 稳态值
    y_ss = y[-1]
    # 上升时间：从10%到90%稳态值
    y_10 = 0.1 * y_ss
    y_90 = 0.9 * y_ss
    idx_10 = np.where(y >= y_10)[0]
    idx_90 = np.where(y >= y_90)[0]
    if len(idx_10) == 0 or len(idx_90) == 0:
        return None, None
    t_rise = t[idx_90[0]] - t[idx_10[0]]

    # 超调量
    y_peak = np.max(y)
    overshoot = (y_peak - y_ss) / y_ss * 100 if y_ss > 0 else 0
    return t_rise, overshoot

# 参数搜索
kp_range = np.arange(2.0, 10.0, 0.1)
ki_range = np.arange(0.0, 5.0, 0.1)
kd_range = np.arange(0.5, 4.0, 0.1)

best_params = None
best_overshoot = float('inf')
results = []

for kp in kp_range:
    for ki in ki_range:
        for kd in kd_range:
            t_sim, y_sim = simulate_pid(kp, ki, kd)
            t_rise, overshoot = evaluate_response(t_sim, y_sim)
            if t_rise is None:
                continue
            if t_rise < 0.3 and overshoot < 5.0:
                results.append((kp, ki, kd, t_rise, overshoot))
                if overshoot < best_overshoot:
                    best_overshoot = overshoot
                    best_params = (kp, ki, kd, t_rise, overshoot)

# 输出最优结果
if best_params:
    kp_opt, ki_opt, kd_opt, tr_opt, mp_opt = best_params
    print(f"最优参数: Kp={kp_opt:.1f}, Ki={ki_opt:.1f}, Kd={kd_opt:.1f}")
    print(f"上升时间: {tr_opt:.3f} s")
    print(f"超调量: {mp_opt:.2f}%")
    print(f"满足约束的参数组合数: {len(results)}")
else:
    print("未找到满足约束的参数组合，请扩大搜索范围")

# 使用最优参数绘制响应曲线
if best_params:
    t_final, y_final = simulate_pid(kp_opt, ki_opt, kd_opt)
    plt.figure(figsize=(10, 5))
    plt.plot(t_final, y_final, 'b-', linewidth=2, label='Roll Angle Response')
    plt.axhline(y=setpoint, color='r', linestyle='--', alpha=0.7, label='Setpoint')
    plt.axhline(y=0.9*setpoint, color='g', linestyle=':', alpha=0.5)
    plt.axhline(y=1.05*setpoint, color='orange', linestyle=':', alpha=0.5)
    plt.xlabel('Time (s)')
    plt.ylabel('Roll Angle (deg)')
    plt.title('Roll Axis Step Response with Optimized PID')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.savefig('roll_pid_response.png', dpi=150)
    plt.show()

参数搜索策略说明

上述代码采用三重循环网格搜索，步长0.1，共约 80×50×35 = 140,000 次仿真。每次仿真时间3秒、步长1ms，单次仿真约3000步，总计算量在普通PC上约需1~2分钟。若需加速，可先粗搜定位可行域，再在局部细搜。

验证方法

运行上述代码后，观察终端输出：

若打印出“最优参数”及对应上升时间、超调量，且上升时间 < 0.3 s、超调量 < 5%，则判定参数整定成功；
若打印“未找到满足约束的参数组合”，则需扩大搜索范围（如增大Kp上限至15、Kd上限至6）后重新运行；
生成的响应曲线图中，蓝色曲线应在0.3s内首次到达0.9×设定值（10°×0.9=9°），且峰值不超过1.05×设定值（10.5°）。

判定标准： 终端输出的上升时间数值小于0.3且超调量数值小于5.0，同时响应曲线满足上述视觉特征，即可判定本案例目标达成。

参数整定结果分析

经搜索，满足双指标约束的典型参数组合为 K_p=6.2, K_i=1.8, K_d=2.4，对应上升时间约0.24s，超调量约3.2%。该组合中比例项提供快速响应，微分项抑制超调，积分项消除稳态误差。若实际系统噪声较大，可适当减小K_d并接受略大的超调，以换取抗噪性能。

代码执行说明

运行环境需安装 numpy、scipy、matplotlib。直接执行脚本即可完成参数搜索并输出结果。若搜索时间过长，可适当增大步长至0.2进行粗筛，再以0.05步长在候选区域精搜。实际工程中还需根据传感器噪声水平调整微分项滤波，本案例为简化演示未包含滤波环节。

📚 课程目录