2. 基于Ziegler-Nichols规则的自动舵PID参数整定与阶跃响应验证

被控对象与整定目标

本例采用一阶Nomoto模型描述航向对舵角的响应：T·ṙ + r = K·δ，其中 r 为艏摇角速度，δ 为舵角，K 为增益，T 为时间常数。取教学示例参数 K = 0.18 1/s、T = 12 s，并加入舵角限幅 ±35°。控制目标是航向阶跃 20° 时：超调量 ≤ 10%、上升时间（10%→90%）≤ 25 s、稳态误差 ≤ 0.5°。

Ziegler-Nichols 临界比例度法的核心是：先关闭积分与微分，仅保留比例环节，逐步增大 Kp 直到系统出现等幅振荡，记录此时的临界增益 Ku 与临界振荡周期 Tu，再按经验公式计算 PID 参数。

临界点搜索与参数计算

对一阶Nomoto模型加纯比例控制，闭环特征方程为 T·s² + s + K·Kp = 0。令 s = jω 代入，实部与虚部分别为零，可得临界条件：ωu = 1/T，Ku = 1/(K·T)，Tu = 2π/ωu = 2πT。代入示例参数：

• Ku = 1/(0.18 × 12) ≈ 0.463
• Tu = 2π × 12 ≈ 75.4 s

按 Z-N 临界比例度法 PID 公式：Kp = 0.6·Ku、Ti = 0.5·Tu、Td = 0.125·Tu，得到 Kp ≈ 0.278、Ti ≈ 37.7 s、Td ≈ 9.4 s，对应 Ki = Kp/Ti ≈ 0.00737、Kd = Kp·Td ≈ 2.62。

离散仿真实现

下面给出完整可运行的 Python 仿真：用欧拉法积分 Nomoto 模型，PID 采用位置式离散形式，并对舵角做 ±35° 限幅。代码同时完成临界点扫描与整定后阶跃响应验证。

import numpy as np

K, T = 0.18, 12.0          # Nomoto 参数
dt = 0.05                  # 仿真步长 (s)
T_end = 300.0
N = int(T_end / dt)
t = np.arange(N) * dt
psi_ref = np.deg2rad(20.0) # 目标航向 20°
delta_max = np.deg2rad(35.0)

def simulate(kp, ki, kd, psi_ref=psi_ref):
    psi, r = 0.0, 0.0
    integ, prev_e = 0.0, 0.0
    hist = np.zeros(N)
    for k in range(N):
        e = psi_ref - psi
        integ += e * dt
        deriv = (e - prev_e) / dt
        delta = kp * e + ki * integ + kd * deriv
        delta = np.clip(delta, -delta_max, delta_max)
        r_dot = (K * delta - r) / T
        r += r_dot * dt
        psi += r * dt
        prev_e = e
        hist[k] = psi
    return hist

def metrics(hist):
    target = psi_ref
    peak = hist.max()
    overshoot = max(0.0, (peak - target) / target * 100)
    # 上升时间 10% -> 90%
    i10 = np.argmax(hist >= 0.1 * target)
    i90 = np.argmax(hist >= 0.9 * target)
    tr = t[i90] - t[i10]
    ess = abs(hist[-1] - target)
    return overshoot, tr, np.rad2deg(ess)

# 1) 临界点扫描：仅比例，找等幅振荡
for kp in np.arange(0.40, 0.55, 0.01):
    h = simulate(kp, 0.0, 0.0)
    amp = h[-50:].max() - h[-50:].min()
    print(f"Kp={kp:.2f}  末段振幅={np.rad2deg(amp):.3f}°")

# 2) 按 Z-N 公式整定
Ku, Tu = 0.463, 2 * np.pi * T
kp = 0.6 * Ku
ti = 0.5 * Tu
td = 0.125 * Tu
ki, kd = kp / ti, kp * td
print(f"Kp={kp:.4f}, Ki={ki:.5f}, Kd={kd:.4f}")

# 3) 阶跃响应验证
h = simulate(kp, ki, kd)
os_, tr_, ess_ = metrics(h)
print(f"超调={os_:.2f}%  上升时间={tr_:.2f}s  稳态误差={ess_:.3f}°")

结果判读与调参方向

按上述参数运行后，典型现象是：超调量偏大（Z-N 公式本身偏激进，常给出 20%~40% 超调），上升时间较快，稳态误差接近零。若超调超过 10% 指标，可按以下方向微调：

一个实用的折中做法是将 Kp 取为 0.45·Ku、Ti 取 0.6·Tu，在保持上升时间满足要求的同时把超调压到 10% 以内。最终是否达标，以仿真输出的三项指标与判定条件逐项比对为准，而不是凭经验公式直接下结论。