用Python“复现”经典控制理论：从传递函数到PID调参的实战指南

如果你正在学习《自动控制原理》，无论是卢京潮老师的课程还是其他经典教材，大概率会经历这样一个阶段：对着满屏的拉普拉斯变换、根轨迹和伯德图感到抽象又困惑。公式推导很严谨，物理概念也清晰，但“这个二阶系统超调量到底是多少？”“这个PID参数调了之后响应曲线会变成什么样？”——这些问题单靠纸笔计算和想象，总感觉隔着一层纱。

作为一名从学生时代过来，也曾在工业现场调试过无数回路的工程师，我深切体会到，将理论“可视化”和“可操作化”是打通任督二脉的关键。而Python，正是我们手中那把绝佳的“手术刀”。它不仅能帮你验证课本上的结论，更能让你亲手“摆弄”系统，观察参数变化的每一个细微影响，把抽象的控制律变成屏幕上直观的曲线。这篇文章，我们就抛开纯理论的复述，直接进入Jupyter Notebook，用代码一行行搭建起从传递函数到控制器设计的完整仿真环境，让你获得一种“亲手掌控系统”的实践感。

1. 环境搭建与核心工具库：打造你的控制仿真实验室

工欲善其事，必先利其器。在开始“控制”之前，我们先得布置好“实验室”。对于控制系统的仿真与分析，Python生态中有几个库堪称神器，它们各有侧重，组合使用能覆盖从建模、分析到控制器设计的全流程。

首先，确保你有一个可运行的Python环境（推荐3.8及以上版本）。使用Anaconda进行环境管理是个好习惯，可以避免库版本冲突。我们核心依赖以下库：

• NumPy & SciPy: 科学计算的基石。NumPy提供高效的数组操作，SciPy则包含了丰富的科学计算模块，其scipy.signal模块是处理线性时不变（LTI）系统的核心，内置了传递函数、状态空间模型的表示以及时域、频域响应的计算函数。
• Matplotlib: 数据可视化的绝对主力。控制系统的性能指标，如阶跃响应曲线、伯德图、奈奎斯特图、根轨迹，最终都要靠它来绘制。
• Control Systems Library (python-control): 这是控制领域的专业库。虽然scipy.signal功能强大，但control库的API设计更贴近控制工程师的思维习惯，尤其在绘制根轨迹、奈奎斯特图、进行系统互联（串联、并联、反馈）等方面更加直观便捷。
• SymPy (可选但推荐): 符号计算库。当你需要推导传递函数、进行拉普拉斯反变换等符号运算时，它能大显身手，帮你验证手工计算的结果。

安装命令非常简单：

pip install numpy scipy matplotlib
pip install control  # 或者使用 conda install -c conda-forge control
pip install sympy

提示：control库在某些系统上可能需要额外安装依赖。如果在导入时遇到问题，可以查阅其官方文档，通常安装 slycot 或 matplotlib 的特定后端可以解决。

下面，我们快速检验一下环境，并创建一个简单的传递函数模型。假设我们有一个经典的二阶系统传递函数：$G(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}$，其中 $\omega_n = 5 , \text{rad/s}$, $\zeta = 0.6$。

import numpy as np
import matplotlib.pyplot as plt
import control as ct

# 定义系统参数
omega_n = 5.0   # 自然频率
zeta = 0.6      # 阻尼比

# 使用 control 库创建传递函数
# 分子多项式系数 (按s的降幂排列): [omega_n**2] 即 [25]
# 分母多项式系数: [1, 2*zeta*omega_n, omega_n**2] 即 [1, 6, 25]
num = [omega_n**2]
den = [1, 2*zeta*omega_n, omega_n**2]
sys_tf = ct.TransferFunction(num, den)

print("系统传递函数 G(s) = ")
print(sys_tf)

# 计算并绘制阶跃响应
t, y = ct.step_response(sys_tf)  # 返回时间序列和输出序列
plt.figure(figsize=(10, 6))
plt.plot(t, y, linewidth=2)
plt.grid(True, which='both', linestyle='--', alpha=0.7)
plt.xlabel('时间 (秒)')
plt.ylabel('输出 y(t)')
plt.title(f'二阶系统阶跃响应 ($\\omega_n$={omega_n}, $\\zeta$={zeta})')
plt.show()

运行这段代码，你将立刻看到一条典型的欠阻尼二阶系统阶跃响应曲线。这就是我们仿真实验室的“第一道光”。接下来，我们将利用这个基础，深入更多经典案例。

2. 从理论到代码：建立系统数学模型

控制系统的分析始于模型。教材中通常会介绍微分方程、传递函数、结构图（方框图）和状态空间方程等多种模型。我们用Python可以轻松地在这些表示形式之间转换和验证。

2.1 传递函数与系统互联

传递函数是最常用的输入输出模型。我们不仅可以直接定义，还可以通过系统互联（串联、并联、反馈）来构建复杂系统。例如，考虑一个前向通道为 $G(s) = \frac{1}{s(s+2)}$，采用单位负反馈的系统。

# 定义前向通道传递函数 G(s)
G = ct.TransferFunction([1], [1, 2, 0])  # 1/(s^2 + 2s) -> 1/(s(s+2))
print("前向通道 G(s):", G)

# 构建单位负反馈闭环系统
sys_cl = ct.feedback(G, 1)  # 第二个参数为1，表示反馈通道传递函数 H(s)=1
print("\n闭环系统传递函数 Φ(s):", sys_cl)

# 比较开环与闭环系统的阶跃响应
t, y_open = ct.step_response(G, T=np.linspace(0, 10, 1000))
t, y_closed = ct.step_response(sys_cl, T=np.linspace(0, 10, 1000))

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(t, y_open, 'r--', linewidth=2, label='开环 G(s)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('时间 (秒)')
plt.ylabel('输出')
plt.title('开环系统阶跃响应')
plt.legend()

plt.subplot(1, 2, 2)
plt.plot(t, y_closed, 'b-', linewidth=2, label='闭环 Φ(s)')
plt.grid(True, linestyle='--', alpha=0.7)
plt.xlabel('时间 (秒)')
plt.ylabel('输出')
plt.title('单位负反馈闭环系统阶跃响应')
plt.legend()
plt.tight_layout()
plt.show()

通过对比，你可以直观地看到反馈如何改变了系统的动态特性（本例中，开环系统相当于一个积分环节，输出会不断增长，而闭环系统则稳定在一个常值）。

2.2 结构图化简与梅森公式验证

对于复杂的多回路系统，手动化简结构图或应用梅森增益公式容易出错。我们可以用control库的series, parallel, feedback函数一步步化简，或者直接定义所有环节和求和点，让库函数自动计算总传递函数。这本质上是在代码中“搭建”结构图。

假设一个稍复杂的系统：$G1(s) = \frac{1}{s+1}$, $G2(s) = \frac{2}{s+3}$， $H(s) = 0.5$， 构成一个前向通道为G1和G2串联，并受到H负反馈的系统。

G1 = ct.TransferFunction([1], [1, 1])
G2 = ct.TransferFunction([2], [1, 3])
H = ct.TransferFunction([0.5], [1])

# 方法1：分步计算
G_forward = ct.series(G1, G2)          # G1 * G2
sys_cl_complex = ct.feedback(G_forward, H) # (G1*G2) / (1 + G1*G2*H)
print("分步计算得到的闭环传递函数:", sys_cl_complex)

# 方法2：验证梅森公式
# 前向通路增益 P1 = G1*G2
# 回路增益 L1 = -G1*G2*H
# 特征式 Δ = 1 - L1 = 1 + G1*G2*H
# 传递函数 Φ = P1 / Δ = (G1*G2) / (1 + G1*G2*H)
# 可见与方法1结果一致。

通过代码验证，可以加深对结构图化简规则和梅森公式的理解，确保理论推导的正确性。

2.3 状态空间模型

现代控制理论基于状态空间模型。control和scipy.signal都支持传递函数与状态空间模型之间的转换。这对于分析系统能控性、能观性以及设计状态反馈控制器至关重要。

# 将之前的二阶系统传递函数转换为状态空间形式
sys_ss = ct.tf2ss(sys_tf)  # 转换为状态空间对象
print("状态空间表示 (A, B, C, D):")
print("A =\n", sys_ss.A)
print("B =\n", sys_ss.B)
print("C =\n", sys_ss.C)
print("D =\n", sys_ss.D)

# 也可以从状态空间方程直接创建
# 例如：一个简单的双积分器系统 x1' = x2, x2' = u, y = x1
A = [[0, 1],
     [0, 0]]
B = [[0],
     [1]]
C = [[1, 0]]
D = [[0]]
sys_integrator = ct.StateSpace(A, B, C, D)
print("\n双积分器系统的传递函数:", ct.ss2tf(sys_integrator))

掌握模型间的转换，意味着你可以在时域和复频域之间自由切换，选择最合适的工具进行分析。

3. 时域分析实战：性能指标与系统响应

时域分析直接反映了系统对输入信号的跟踪能力。我们最关心的是稳定性、快速性（调节时间、上升时间）和平稳性（超调量）。

3.1 阶跃响应与性能指标提取

我们可以仿真阶跃响应，并编写函数自动计算关键指标。

def step_response_metrics(sys, T_end=10, plot=True):
    """计算并显示系统阶跃响应的关键性能指标"""
    t, y = ct.step_response(sys, T=np.linspace(0, T_end, 5000))
    y_final = y[-1]

    # 1. 稳态值
    y_ss = y_final

    # 2. 上升时间 (从10%到90%)
    idx_10 = np.where(y >= 0.1 * y_ss)[0][0]
    idx_90 = np.where(y >= 0.9 * y_ss)[0][0]
    t_r = t[idx_90] - t[idx_10]

    # 3. 峰值时间与超调量
    y_max = np.max(y)
    t_p = t[np.argmax(y)]
    overshoot = (y_max - y_ss) / y_ss * 100 if y_ss != 0 else 0

    # 4. 调节时间 (进入并保持在±2%误差带内的时间)
    # 寻找最后一个超出误差带的时间点
    err_band = 0.02
    within_band = np.abs(y - y_ss) <= err_band * np.abs(y_ss)
    # 从后往前找，第一个为False的索引之后的时间都是满足条件的
    settling_indices = np.where(~within_band)[0]
    t_s = t[settling_indices[-1] + 1] if len(settling_indices) > 0 else 0

    metrics = {
        '稳态值': y_ss,
        '上升时间 (10%-90%)': t_r,
        '峰值时间': t_p,
        '最大超调量 (%)': overshoot,
        '调节时间 (±2%)': t_s
    }

    if plot:
        plt.figure(figsize=(10, 6))
        plt.plot(t, y, 'b-', linewidth=2, label='阶跃响应')
        plt.axhline(y=y_ss, color='k', linestyle='--', alpha=0.5, label='稳态值')
        plt.axhline(y=y_ss*(1+err_band), color='r', linestyle=':', alpha=0.7, label='±2%误差带')
        plt.axhline(y=y_ss*(1-err_band), color='r', linestyle=':', alpha=0.7)
        plt.axvline(x=t_p, color='g', linestyle='--', alpha=0.7, label=f'峰值时间 t_p={t_p:.2f}s')
        plt.fill_between(t, y_ss*(1-err_band), y_ss*(1+err_band), color='gray', alpha=0.1)
        plt.grid(True, linestyle='--', alpha=0.7)
        plt.xlabel('时间 (秒)')
        plt.ylabel('输出 y(t)')
        plt.title('系统阶跃响应及性能指标')
        plt.legend()
        plt.show()

    return metrics, t, y

# 测试我们的函数
metrics, _, _ = step_response_metrics(sys_cl)  # 使用之前定义的闭环系统
print("\n系统性能指标：")
for key, value in metrics.items():
    print(f"  {key}: {value:.4f}")

这个自定义函数不仅绘制了响应曲线，还自动标注了关键指标，让你对系统性能一目了然。你可以修改系统参数（如阻尼比ζ），观察这些指标如何变化，从而深刻理解公式 $t_r \approx \frac{1.8}{\omega_n}$, $M_p = e^{-\pi\zeta/\sqrt{1-\zeta^2}} \times 100%$, $t_s \approx \frac{4.6}{\zeta\omega_n}$ (2%准则) 的物理意义。

3.2 系统稳定性判定：劳斯-赫尔维茨判据

虽然对于低阶系统，直接求特征根更直观，但劳斯表对于高阶系统或含参数的系统稳定性分析非常有用。我们可以用control库的routh函数（或scipy.signal的routh）来判定。

# 定义一个三阶系统：s^3 + 3s^2 + 3s + K，分析K对稳定性的影响
def check_stability_by_routh(K):
    """使用劳斯判据判断系统稳定性"""
    den = [1, 3, 3, K]  # 特征多项式系数
    # 注意：control.routh 返回劳斯表，但需要手动判断第一列符号变化
    # 这里我们使用一个更直接的方法：求根
    roots = np.roots(den)
    # 判断是否有实部大于等于0的根（临界稳定算不稳定）
    unstable = any(np.real(root) >= 0 for root in roots)
    return not unstable, roots

# 测试不同的K值
K_values = [1, 5, 10, 15]
stability_table = []
for K in K_values:
    stable, roots = check_stability_by_routh(K)
    stability_table.append([K, stable, roots])

# 用表格展示结果
import pandas as pd
df = pd.DataFrame(stability_table, columns=['K值', '是否稳定', '特征根'])
print(df.to_string(index=False))

通过改变K值，你会发现当K>9时，系统变得不稳定（有正实部根）。这比单纯记忆劳斯表第一列全为正的判据要生动得多。

4. 频域分析与综合：伯德图、奈奎斯特图与校正设计

频域分析是经典控制的精华，它通过开环频率特性来预测闭环系统的性能（稳定性、稳态误差、动态响应）。

4.1 绘制伯德图与奈奎斯特图

control库使绘制这些图变得异常简单。

# 继续使用之前的闭环系统开环传递函数 G(s) = 1/(s(s+2))
plt.figure(figsize=(14, 5))

# 伯德图
plt.subplot(1, 2, 1)
mag, phase, omega = ct.bode_plot(G, dB=True, Hz=False, margins=True, omega_limits=(0.1, 100))
plt.grid(True, which='both', linestyle='--', alpha=0.7)
plt.title('开环系统伯德图')

# 奈奎斯特图
plt.subplot(1, 2, 2)
ct.nyquist_plot(G, omega_limits=(0.1, 100), arrows=5)
plt.grid(True, linestyle='--', alpha=0.7)
plt.title('开环系统奈奎斯特图')
plt.tight_layout()
plt.show()

# 从伯德图中读取幅值裕度和相位裕度
gm, pm, wpc, wgc = ct.margin(G)
print(f"幅值裕度 Gm = {gm:.2f} (即 {20*np.log10(gm):.2f} dB) 在频率 ω = {wpc:.2f} rad/s")
print(f"相位裕度 Pm = {pm:.2f}° 在频率 ω = {wgc:.2f} rad/s")

伯德图上的幅值裕度和相位裕度是衡量系统相对稳定性的关键指标。通过代码，你可以精确地获取它们的数值，并理解“裕度越大，系统相对稳定性越好”的含义。

4.2 基于频域的串联校正设计

当系统性能不满足要求时，需要设计校正环节。以串联超前校正为例，其传递函数为 $G_c(s) = K_c \frac{Ts+1}{\alpha Ts+1}, , \alpha < 1$。它的核心思想是利用其提供的正相位补偿，提高系统的相位裕度，从而改善动态响应。

假设我们的原系统 $G(s) = \frac{1}{s(s+1)}$，要求静态速度误差系数 $K_v \geq 20$，相位裕度 $\gamma \geq 50°$。

# 1. 原系统分析
G_original = ct.TransferFunction([1], [1, 1, 0])  # 1/(s(s+1))
print("原系统开环传递函数:", G_original)

# 计算原系统在Kv=20时的增益K
# Kv = lim_{s->0} s * K * G(s) = K * 1 = K， 所以 K = 20
K = 20
G_original_scaled = ct.TransferFunction([K], [1, 1, 0])

# 计算此时的相位裕度
_, pm_orig, _, wgc_orig = ct.margin(G_original_scaled)
print(f"增益调整后 (K={K}) 的原系统相位裕度: {pm_orig:.1f}° (目标≥50°)")

# 2. 设计超前校正网络
# 需要补偿的相位增量 φm = γ_target - γ_original + (5°~12°补偿余量)
gamma_target = 50
compensation = 10  # 补偿余量
phi_m = np.deg2rad(gamma_target - pm_orig + compensation)

# 计算超前网络的参数 α
alpha = (1 - np.sin(phi_m)) / (1 + np.sin(phi_m))
print(f"所需最大超前相位 φm = {np.rad2deg(phi_m):.1f}°")
print(f"计算得到的 α = {alpha:.3f}")

# 确定校正网络的两个转折频率：将最大超前相位处的频率 ωm 置于新的剪切频率处
# 我们需要找到一个新的剪切频率 ωgc_new，使得在此频率下，原系统的幅值为 -10log10(1/α) dB
# 即 |G(jωgc_new)| = 1/sqrt(α)  (线性尺度)
# 这里我们通过迭代或查找伯德图数据来近似求解
mag, phase, omega = ct.bode_plot(G_original_scaled, plot=False, omega=np.logspace(-1, 2, 1000))
mag_linear = 10**(mag/20)  # 转换为线性幅值
# 寻找幅值最接近 1/np.sqrt(alpha) 的频率点
target_mag = 1/np.sqrt(alpha)
idx = np.argmin(np.abs(mag_linear - target_mag))
omega_m = omega[idx]
print(f"预估的 ωm (新剪切频率附近) = {omega_m:.2f} rad/s")

# 计算时间常数 T
T = 1 / (omega_m * np.sqrt(alpha))
print(f"计算得到的时间常数 T = {T:.3f}")

# 构建超前校正网络传递函数
Gc_lead = ct.TransferFunction([K * T, K], [alpha * T, 1])  # K * (Ts+1)/(αTs+1)
print(f"超前校正网络传递函数: {Gc_lead}")

# 3. 校正后系统分析
G_corrected = ct.series(Gc_lead, G_original)  # Gc(s) * G(s)
print("\n校正后开环传递函数:", G_corrected)

# 绘制校正前后伯德图对比
plt.figure(figsize=(12, 8))
mag_orig, phase_orig, omega_orig = ct.bode_plot(G_original_scaled, plot=False, omega=np.logspace(-1, 2, 1000))
mag_corr, phase_corr, omega_corr = ct.bode_plot(G_corrected, plot=False, omega=np.logspace(-1, 2, 1000))

plt.subplot(2, 1, 1)
plt.semilogx(omega_orig, 20*np.log10(mag_orig), 'b--', label='原系统 (调整增益后)', linewidth=2)
plt.semilogx(omega_corr, 20*np.log10(mag_corr), 'r-', label='校正后系统', linewidth=2)
plt.grid(True, which='both', linestyle='--', alpha=0.7)
plt.ylabel('幅值 (dB)')
plt.title('超前校正前后伯德图对比')
plt.legend()

plt.subplot(2, 1, 2)
plt.semilogx(omega_orig, phase_orig, 'b--', label='原系统', linewidth=2)
plt.semilogx(omega_corr, phase_corr, 'r-', label='校正后系统', linewidth=2)
plt.grid(True, which='both', linestyle='--', alpha=0.7)
plt.ylabel('相位 (度)')
plt.xlabel('频率 (rad/s)')
plt.legend()
plt.tight_layout()
plt.show()

# 验证校正后系统性能
gm_corr, pm_corr, wpc_corr, wgc_corr = ct.margin(G_corrected)
print(f"\n校正后系统性能验证：")
print(f"  相位裕度 Pm = {pm_corr:.1f}° (目标≥50°)")
# 验证静态速度误差系数
Kv_corrected = ct.dcgain(ct.TransferFunction([1, 0], [1]) * G_corrected)  # s * G(s) |_{s=0}
print(f"  静态速度误差系数 Kv = {Kv_corrected:.1f} (目标≥20)")

这段代码完整演示了基于频率响应法的超前校正设计流程：从性能指标分析、计算校正网络参数，到验证校正效果。你可以调整目标相位裕度或原系统，观察校正网络参数如何变化，并最终通过伯德图和阶跃响应验证设计是否成功。这种“设计-验证”的闭环体验，是单纯看书做题难以获得的。

5. PID控制器整定与仿真：当理论遇上实践

PID控制器是工业应用的绝对主力。理论课程会讲齐格勒-尼古拉斯（Z-N）整定法等，但实际调参往往更依赖经验和试凑。我们可以用仿真来模拟这个过程，理解P、I、D三个参数各自的作用。

5.1 模拟一个被控对象并手动整定PID

假设我们控制一个带延迟的一阶惯性环节：$G_p(s) = \frac{K e^{-\theta s}}{\tau s + 1}$， 其中 $K=2, \tau=5, \theta=1$。我们使用一个简单的离散PID算法进行仿真。

import matplotlib.pyplot as plt
import numpy as np

def simulate_pid(Kp, Ki, Kd, setpoint=1.0, T=50, dt=0.1):
    """
    离散PID控制仿真
    被控对象: 一阶惯性加纯延迟 (用一阶帕德近似表示延迟)
    """
    # 连续对象模型: G(s) = 2 * exp(-1*s) / (5s + 1)
    # 用一阶帕德近似 exp(-theta*s) ≈ (1 - theta*s/2) / (1 + theta*s/2)
    # 因此近似对象传递函数为: G_app(s) = 2 * (1 - 0.5*s) / (5s+1)(1+0.5s)
    # 我们将其转换为状态空间或直接进行离散化仿真
    # 这里为了简化，用差分方程模拟一阶惯性，并加上一个步长的延迟
    tau = 5.0
    K = 2.0
    delay_steps = int(1.0 / dt)  # 1秒延迟对应的步数

    n_steps = int(T / dt)
    time = np.arange(0, T, dt)
    y = np.zeros(n_steps)  # 过程输出
    u = np.zeros(n_steps)  # 控制器输出
    e = np.zeros(n_steps)  # 误差
    e_int = 0.0  # 误差积分
    e_prev = 0.0  # 上一时刻误差 (用于微分)

    # 初始化（假设初始状态为0）
    y_hist = np.zeros(delay_steps)  # 用于存储延迟的输出历史

    for i in range(1, n_steps):
        # 计算误差
        e[i] = setpoint - y[i-1]

        # PID计算 (位置式)
        e_int += e[i] * dt
        e_deriv = (e[i] - e_prev) / dt if dt > 0 else 0

        u[i] = Kp * e[i] + Ki * e_int + Kd * e_deriv
        # 简单的输出限幅
        u[i] = np.clip(u[i], 0, 5)

        # 被控对象动态 (一阶惯性离散化)
        # dy/dt = (K*u - y) / tau
        y_dot = (K * u[i] - y[i-1]) / tau
        y[i] = y[i-1] + y_dot * dt

        # 加入纯延迟：当前输出实际上是 delay_steps 步之前的输入作用的结果
        if i >= delay_steps:
            # 这里简化处理，将延迟作用在输出上（更准确的应作用在输入到对象）
            # 使用历史缓冲区模拟
            y[i] = y_hist[0]  # 取最早的历史值作为当前输出
            y_hist = np.roll(y_hist, -1)  # 历史缓冲区向前移动
            y_hist[-1] = y[i]  # 将当前计算值放入缓冲区末尾（用于未来输出）
        else:
            # 初始阶段，延迟未满，直接使用当前计算值
            y_hist[i-1] = y[i]

        e_prev = e[i]

    # 绘制结果
    plt.figure(figsize=(12, 8))
    plt.subplot(2, 1, 1)
    plt.plot(time, [setpoint]*len(time), 'k--', label='设定值')
    plt.plot(time, y, 'b-', linewidth=2, label='过程输出 (y)')
    plt.grid(True, linestyle='--', alpha=0.7)
    plt.ylabel('输出')
    plt.title(f'PID控制响应 (Kp={Kp}, Ki={Ki}, Kd={Kd})')
    plt.legend()
    plt.xlim([0, T])

    plt.subplot(2, 1, 2)
    plt.plot(time, u, 'r-', linewidth=2, label='控制量 (u)')
    plt.grid(True, linestyle='--', alpha=0.7)
    plt.xlabel('时间 (秒)')
    plt.ylabel('控制量')
    plt.legend()
    plt.xlim([0, T])
    plt.tight_layout()
    plt.show()

    # 计算性能指标
    y_steady = y[-int(5/dt):].mean()  # 最后5秒的平均作为稳态值
    overshoot = (np.max(y) - setpoint) / setpoint * 100 if setpoint != 0 else 0
    # 粗略计算调节时间 (进入±2%误差带)
    err_band = 0.02 * setpoint
    idx_settled = np.where(np.abs(y - setpoint) <= err_band)[0]
    t_settle = time[idx_settled[0]] if len(idx_settled) > 0 else T

    print(f"性能指标:")
    print(f"  稳态值: {y_steady:.3f} (设定值: {setpoint})")
    print(f"  超调量: {overshoot:.1f}%")
    print(f"  调节时间(近似): {t_settle:.1f} s")

# 尝试不同的PID参数
print("=== 纯比例控制 (P) ===")
simulate_pid(Kp=1.0, Ki=0.0, Kd=0.0)

print("\n=== 比例积分控制 (PI) ===")
simulate_pid(Kp=1.2, Ki=0.2, Kd=0.0)

print("\n=== 完整的PID控制 ===")
simulate_pid(Kp=1.5, Ki=0.25, Kd=0.8)

运行这段代码，你会看到三组不同参数下的控制效果。纯比例控制会有稳态误差；加入积分作用可以消除稳态误差，但可能引起超调和振荡；再加入微分作用，可以预测误差变化趋势，抑制超调，提高响应速度。通过反复调整参数并观察响应曲线，你能直观地感受到每个参数是如何影响系统动态的，这正是PID整定的核心。

5.2 利用仿真优化PID参数

手动试凑效率低。我们可以结合一些启发式规则或简单的优化算法来辅助整定。例如，先使用齐格勒-尼古拉斯临界比例度法确定一组初始参数，再进行微调。

# 假设我们通过实验或仿真，得到了对象的临界增益Ku和临界振荡周期Tu
# 例如，对于上述对象，通过仿真发现当Kp≈3.2时，系统出现等幅振荡（临界稳定），振荡周期Tu≈8s
Ku = 3.2
Tu = 8.0

# Z-N整定公式 (PID型)
Kp_zn = 0.6 * Ku
Ti_zn = 0.5 * Tu
Td_zn = 0.125 * Tu
Ki_zn = Kp_zn / Ti_zn
Kd_zn = Kp_zn * Td_zn

print(f"基于Z-N整定法的PID参数：")
print(f"  Kp = {Kp_zn:.3f}")
print(f"  Ki = {Ki_zn:.3f} (Ti={Ti_zn:.2f}s)")
print(f"  Kd = {Kd_zn:.3f} (Td={Td_zn:.2f}s)")

print("\n=== Z-N整定PID控制效果 ===")
simulate_pid(Kp=Kp_zn, Ki=Ki_zn, Kd=Kd_zn)

将Z-N法得到的参数代入仿真，你可能会得到一个可用的控制器，但性能未必最优。这时，可以在其附近进行参数扫描，寻找更优解。这种基于仿真的“虚拟调试”，成本极低，却能极大加深你对PID控制器和被控对象之间相互作用的理解。

从搭建仿真环境，到建立系统模型，再到时域频域分析和控制器设计，最后进行PID整定，我们完成了一个完整的控制理论学习与实践的闭环。代码不仅仅是验证理论的工具，它更是一个沙盘，让你可以大胆尝试各种“如果”：如果阻尼比变小会怎样？如果增益变大系统会失稳吗？这个校正网络真的有效吗？这些问题的答案，不再停留在公式和想象中，而是变成了一条条你可以亲眼观察、亲手测量的曲线。