Python实战：用control库5分钟搞定Bode图与Nyquist图绘制（附完整代码）

如果你正在学习自动控制原理，或者工作中需要分析线性系统的频率响应，那么Bode图和Nyquist图这两个工具你一定不陌生。它们就像是工程师的“听诊器”，能让我们“听”到系统在不同频率下的行为特征。但说实话，传统的手工绘制方法既繁琐又容易出错，特别是面对复杂的传递函数时，计算量让人头疼。

我记得自己刚开始接触频域分析时，光是画一个三阶系统的Nyquist图就要花上大半天，各种角度计算、幅值估算，最后画出来的曲线还常常不够精确。直到后来发现了Python的control库，整个分析过程才变得高效而优雅。现在，我只需要几行代码，就能在几分钟内生成专业级的频率特性图，还能自动计算稳定裕度、带宽等关键指标。

这篇文章就是为你准备的。无论你是自动化专业的学生需要完成课程作业，还是工程师需要快速验证系统设计，我都会带你用最直接的方式掌握Python频域分析的核心技巧。我们不谈复杂的理论推导，只聚焦于如何用代码解决实际问题。我会分享一些我在实际项目中积累的经验，比如如何处理数值计算中的陷阱，如何让图表更清晰美观，以及如何解读这些图形背后的工程意义。

1. 环境搭建与control库快速入门

在开始绘制图表之前，我们需要先搭建好Python环境。我推荐使用Anaconda发行版，它集成了科学计算所需的大部分库，管理起来非常方便。如果你已经安装了Python，也可以通过pip直接安装必要的包。

1.1 安装必要的库

打开你的终端或命令提示符，执行以下命令：

pip install control matplotlib numpy scipy

这里我们主要需要四个库：

• control：控制系统分析与设计的核心库
• matplotlib：绘图库，用于生成Bode图、Nyquist图等
• numpy：数值计算基础库
• scipy：科学计算库，control库有时会依赖它

安装完成后，我们可以通过一个简单的测试来验证是否安装成功：

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

print(f"control库版本: {control.__version__}")
print(f"numpy版本: {np.__version__}")

如果一切正常，你应该能看到各个库的版本号。我目前使用的是control 0.9.4，这个版本已经相当稳定，包含了我们需要的所有频域分析功能。

1.2 创建你的第一个传递函数

在控制理论中，我们通常用传递函数来描述线性时不变系统。control库提供了多种创建传递函数的方式，最常用的是TransferFunction类。

让我们从一个简单的二阶系统开始：

# 创建一个二阶系统：G(s) = 1/(s^2 + 2ζω_n s + ω_n^2)
# 其中阻尼比ζ=0.5，自然频率ω_n=1 rad/s

zeta = 0.5  # 阻尼比
omega_n = 1.0  # 自然频率

# 传递函数的分母多项式系数：s^2 + 2*zeta*omega_n*s + omega_n^2
den = [1, 2*zeta*omega_n, omega_n**2]
num = [1]  # 分子为1

sys = control.TransferFunction(num, den)
print(f"系统传递函数: {sys}")

运行这段代码，你会看到输出：

系统传递函数: 
     1
---------------
s^2 + s + 1

提示：在control库中，多项式系数按照降幂排列。比如[1, 2, 1]表示s^2 + 2s + 1。这一点与MATLAB的习惯一致，但如果你之前主要用符号计算库，可能需要适应一下。

1.3 理解频率响应的基本概念

在深入代码之前，我们先快速回顾一下频率响应的核心概念。当你给一个线性系统输入一个正弦信号时，系统的稳态输出也是一个同频率的正弦信号，但幅值和相位会发生变化。这种变化与输入信号的频率有关：

• 幅频特性：输出幅值与输入幅值的比值随频率的变化
• 相频特性：输出相位与输入相位的差值随频率的变化

Bode图就是将这两个特性分别画在对数坐标上，而Nyquist图则是将频率响应画在复平面上。理解这一点很重要，因为后续的所有分析都建立在这个基础之上。

2. Bode图绘制：从基础到高级技巧

Bode图是频域分析中最常用的工具之一，它由两个子图组成：幅频特性图（对数坐标）和相频特性图。control库提供了非常便捷的绘制函数，但要想得到专业级的图表，还需要掌握一些技巧。

2.1 基础Bode图绘制

让我们从最简单的开始。使用前面创建的二阶系统：

# 设置中文字体显示（如果你的系统支持中文）
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 创建图形
plt.figure(figsize=(10, 8))

# 绘制Bode图
mag, phase, omega = control.bode(sys, dB=True, Hz=False, deg=True, plot=True)

# 添加网格和标题
plt.suptitle('二阶系统Bode图 (ζ=0.5, ω_n=1 rad/s)', fontsize=14)
plt.tight_layout()
plt.show()

这段代码会生成一个标准的Bode图。control.bode()函数有几个重要参数需要了解：

参数	说明	默认值
dB	幅值是否用分贝表示	True
Hz	频率单位是否为Hz	False
deg	相位单位是否为度	True
plot	是否自动绘图	True
omega	自定义频率范围	自动生成

注意：默认的频率范围是自动计算的，但对于某些特殊系统（如纯积分或纯微分环节），可能需要手动指定频率范围，否则可能无法正确显示。

2.2 自定义频率范围和分辨率

在实际工程中，我们往往需要关注特定频率范围内的系统行为。比如，在分析伺服系统时，我们可能更关心1-100Hz范围内的特性。control库允许我们自定义频率点：

# 方法1：使用logspace生成对数均匀分布的点
omega_custom = np.logspace(-2, 2, 500)  # 从10^-2到10^2 rad/s，共500个点

# 方法2：使用linspace生成线性均匀分布的点（不推荐，因为Bode图通常用对数坐标）
# omega_custom = np.linspace(0.01, 100, 1000)

# 绘制自定义频率范围的Bode图
plt.figure(figsize=(12, 8))

# 先计算频率响应，不立即绘图
mag, phase, omega = control.bode(sys, omega_custom, dB=True, deg=True, plot=False)

# 手动绘制，以便更精细地控制
plt.subplot(2, 1, 1)
plt.semilogx(omega, 20*np.log10(mag), 'b-', linewidth=2)
plt.grid(True, which='both', linestyle='--', alpha=0.5)
plt.ylabel('幅值 (dB)', fontsize=12)
plt.title('自定义频率范围的Bode图', fontsize=14)

# 添加截止频率标记
mag_db = 20*np.log10(mag)
# 找到幅值穿越0dB的点（截止频率）
zero_crossings = np.where(np.diff(np.sign(mag_db)))[0]
if len(zero_crossings) > 0:
    idx = zero_crossings[0]
    wc = omega[idx]
    plt.axvline(x=wc, color='r', linestyle='--', alpha=0.7)
    plt.text(wc, 5, f'ω_c = {wc:.2f} rad/s', 
             rotation=90, verticalalignment='bottom')

plt.subplot(2, 1, 2)
plt.semilogx(omega, phase, 'g-', linewidth=2)
plt.grid(True, which='both', linestyle='--', alpha=0.5)
plt.xlabel('频率 (rad/s)', fontsize=12)
plt.ylabel('相位 (度)', fontsize=12)

# 添加-180度线（稳定性参考）
plt.axhline(y=-180, color='r', linestyle='--', alpha=0.5)

plt.tight_layout()
plt.show()

这种手动绘制的方式虽然代码量稍多，但给了我们更大的灵活性。比如，我经常需要：

• 标记特定的频率点（如截止频率、谐振频率）
• 添加参考线（如0dB线、-180度线）
• 使用不同的线型和颜色区分多个系统
• 添加图例和注释

2.3 多系统对比分析

在实际工作中，我们经常需要比较不同参数或不同结构系统的频率特性。control库可以轻松实现多系统Bode图的对比：

# 创建不同阻尼比的二阶系统进行对比
damping_ratios = [0.2, 0.5, 0.7, 1.0]  # 不同的阻尼比
omega_n = 2.0  # 固定自然频率

plt.figure(figsize=(12, 10))
colors = ['b', 'g', 'r', 'm']  # 不同颜色
line_styles = ['-', '--', '-.', ':']  # 不同线型

# 幅频特性图
plt.subplot(2, 1, 1)
for i, zeta in enumerate(damping_ratios):
    den = [1, 2*zeta*omega_n, omega_n**2]
    sys_temp = control.TransferFunction([1], den)
    
    mag, phase, omega = control.bode(sys_temp, dB=True, deg=True, plot=False)
    
    plt.semilogx(omega, 20*np.log10(mag), 
                 color=colors[i], 
                 linestyle=line_styles[i],
                 linewidth=2,
                 label=f'ζ={zeta}')

plt.grid(True, which='both', linestyle='--', alpha=0.5)
plt.ylabel('幅值 (dB)', fontsize=12)
plt.title('不同阻尼比二阶系统的Bode图对比', fontsize=14)
plt.legend(loc='best')

# 相频特性图
plt.subplot(2, 1, 2)
for i, zeta in enumerate(damping_ratios):
    den = [1, 2*zeta*omega_n, omega_n**2]
    sys_temp = control.TransferFunction([1], den)
    
    mag, phase, omega = control.bode(sys_temp, dB=True, deg=True, plot=False)
    
    plt.semilogx(omega, phase, 
                 color=colors[i], 
                 linestyle=line_styles[i],
                 linewidth=2,
                 label=f'ζ={zeta}')

plt.grid(True, which='both', linestyle='--', alpha=0.5)
plt.xlabel('频率 (rad/s)', fontsize=12)
plt.ylabel('相位 (度)', fontsize=12)
plt.legend(loc='best')

plt.tight_layout()
plt.show()

通过这样的对比，我们可以直观地看到阻尼比对系统频率特性的影响：

• 阻尼比越小，谐振峰值越大
• 阻尼比影响相位曲线的形状
• 不同阻尼比下，系统的带宽也不同

2.4 自动计算稳定裕度

稳定裕度是衡量系统相对稳定性的重要指标。control库提供了margin()函数，可以自动计算增益裕度、相位裕度以及对应的频率：

# 创建一个三阶系统示例
num = [10]
den = [1, 3, 2, 0]  # s(s+1)(s+2) = s^3 + 3s^2 + 2s
sys_example = control.TransferFunction(num, den)

# 计算稳定裕度
gm, pm, wg, wc = control.margin(sys_example)

print("=== 稳定裕度分析 ===")
print(f"增益裕度 GM = {gm:.4f} (线性值)")
print(f"增益裕度 GM = {20*np.log10(gm):.4f} dB")
print(f"相位裕度 PM = {pm:.4f} 度")
print(f"相位穿越频率 ω_g = {wg:.4f} rad/s")
print(f"截止频率 ω_c = {wc:.4f} rad/s")

# 在Bode图上标注稳定裕度
plt.figure(figsize=(10, 8))
mag, phase, omega = control.bode(sys_example, dB=True, deg=True, plot=True)

# 获取当前坐标轴
ax1, ax2 = plt.gcf().axes

# 在幅频图上标注增益裕度
ax1.axvline(x=wg, color='r', linestyle='--', alpha=0.7)
ax1.axhline(y=20*np.log10(gm), color='r', linestyle='--', alpha=0.7)
ax1.plot(wg, 20*np.log10(gm), 'ro', markersize=8)
ax1.text(wg, 20*np.log10(gm)+2, 
         f'GM={20*np.log10(gm):.1f} dB\nω_g={wg:.2f} rad/s',
         ha='center', fontsize=10)

# 在相频图上标注相位裕度
ax2.axvline(x=wc, color='g', linestyle='--', alpha=0.7)
ax2.axhline(y=-180+pm, color='g', linestyle='--', alpha=0.7)
ax2.plot(wc, -180+pm, 'go', markersize=8)
ax2.text(wc, -180+pm+10, 
         f'PM={pm:.1f}°\nω_c={wc:.2f} rad/s',
         ha='center', fontsize=10)

plt.suptitle('系统Bode图与稳定裕度标注', fontsize=14)
plt.tight_layout()
plt.show()

在实际项目中，我通常会用表格记录不同设计方案的稳定裕度，方便对比：

设计方案	相位裕度 (度)	增益裕度 (dB)	截止频率 (rad/s)	评价
原始设计	15.2	4.8	2.1	裕度不足
超前校正	45.3	12.1	3.5	良好
滞后校正	55.1	8.7	1.2	裕度充足但带宽小
滞后-超前	48.7	10.5	2.8	平衡性好

这样的表格能让设计决策更加直观。通常，工业控制系统要求相位裕度在30-60度之间，增益裕度大于6dB，这样的系统既有足够的稳定性，又有较好的动态响应。

3. Nyquist图绘制与稳定性分析

Nyquist图是另一种重要的频域分析工具，它将频率响应画在复平面上，特别适合用于稳定性分析。与Bode图相比，Nyquist图能更直观地展示系统是否包围(-1, j0)点，这是判断闭环系统稳定性的关键。

3.1 基础Nyquist图绘制

control库提供了nyquist_plot()函数来绘制Nyquist图，但默认设置可能不适合所有情况。让我们从一个简单的例子开始：

# 创建几个不同稳定性的系统进行对比
systems = []

# 稳定系统：G(s) = 1/(s^2 + s + 1)
sys_stable = control.TransferFunction([1], [1, 1, 1])
systems.append(('稳定系统', sys_stable))

# 不稳定系统：G(s) = 1/(s^2 - s - 1) 有一个右半平面极点
sys_unstable = control.TransferFunction([1], [1, -1, -1])
systems.append(('不稳定系统', sys_unstable))

# 临界稳定系统：G(s) = 1/(s^2 + 1) 在虚轴上有极点
sys_marginally = control.TransferFunction([1], [1, 0, 1])
systems.append(('临界稳定系统', sys_marginally))

# 绘制Nyquist图
plt.figure(figsize=(15, 5))

for i, (name, sys) in enumerate(systems):
    plt.subplot(1, 3, i+1)
    
    # 绘制Nyquist图
    control.nyquist_plot(sys, omega_limits=(0.01, 10), 
                         arrow_width=0.02, arrow_color='red')
    
    # 标记(-1, j0)点
    plt.plot(-1, 0, 'ro', markersize=10, markeredgewidth=2, 
             markerfacecolor='none', label='(-1, j0)点')
    
    # 添加单位圆
    theta = np.linspace(0, 2*np.pi, 100)
    plt.plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.3, linewidth=1)
    
    plt.grid(True, alpha=0.5)
    plt.axhline(y=0, color='k', linestyle='-', alpha=0.3, linewidth=0.5)
    plt.axvline(x=0, color='k', linestyle='-', alpha=0.3, linewidth=0.5)
    plt.xlabel('实部', fontsize=11)
    plt.ylabel('虚部', fontsize=11)
    plt.title(f'{name} Nyquist图', fontsize=12)
    plt.legend(loc='best', fontsize=9)
    plt.axis('equal')

plt.tight_layout()
plt.show()

这段代码展示了三个不同稳定性系统的Nyquist图。注意几个关键点：

1. 稳定系统：Nyquist曲线不包围(-1, j0)点
2. 不稳定系统：需要根据开环右半平面极点数判断
3. 临界稳定系统：曲线通过(-1, j0)点

nyquist_plot()函数有几个实用参数：

• omega_limits：指定频率范围，避免数值问题
• arrow_width：箭头宽度，指示频率增加方向
• arrow_color：箭头颜色
• plot：是否立即绘图（默认为True）

3.2 手动实现Nyquist判据

虽然control库没有直接提供Nyquist判据的计算函数，但我们可以基于频率响应数据自己实现。这对于理解判据原理很有帮助：

def nyquist_criterion(sys, omega_range=(0.01, 100), num_points=10000):
    """
    手动实现Nyquist稳定判据
    
    参数:
        sys: 开环传递函数
        omega_range: 频率范围 (rad/s)
        num_points: 频率点数
    
    返回:
        encirclements: 包围(-1, j0)点的圈数
        stability: 闭环稳定性 ('稳定', '不稳定', '临界稳定')
    """
    # 生成频率点（避免ω=0）
    omega = np.logspace(np.log10(omega_range[0]), 
                       np.log10(omega_range[1]), 
                       num_points)
    
    # 计算频率响应
    freqs = 1j * omega
    response = sys(1j * omega)  # 计算频率响应
    
    # 提取实部和虚部
    real_part = np.real(response)
    imag_part = np.imag(response)
    
    # 计算包围(-1, j0)点的圈数
    # 使用角度变化法
    angles = np.angle(response + 1)  # 相对于(-1, j0)点的角度
    angle_change = np.sum(np.diff(np.unwrap(angles)))
    
    # 计算圈数（每2π为一圈）
    encirclements = angle_change / (2 * np.pi)
    
    # 计算开环右半平面极点数
    poles = control.pole(sys)
    P = np.sum(np.real(poles) > 0)
    
    # 计算闭环右半平面极点数
    Z = P - encirclements
    
    # 判断稳定性
    if Z > 0:
        stability = '不稳定'
    elif Z == 0:
        stability = '稳定'
    else:
        stability = '临界稳定'
    
    return {
        'encirclements': encirclements,
        'P': P,
        'Z': Z,
        'stability': stability,
        'real_part': real_part,
        'imag_part': imag_part,
        'omega': omega
    }

# 测试几个系统
test_systems = [
    ('G(s)=1/(s+1)', control.TransferFunction([1], [1, 1])),
    ('G(s)=1/(s^2+s+1)', control.TransferFunction([1], [1, 1, 1])),
    ('G(s)=1/(s-1)', control.TransferFunction([1], [1, -1])),
    ('G(s)=1/(s(s+1))', control.TransferFunction([1], [1, 1, 0])),
]

print("=== Nyquist判据分析结果 ===")
print("-" * 60)
print(f"{'系统':<20} {'P(右半平面极点)':<20} {'包围圈数':<15} {'Z(闭环右极点)':<20} {'稳定性':<10}")
print("-" * 60)

for name, sys in test_systems:
    result = nyquist_criterion(sys)
    print(f"{name:<20} {result['P']:<20.0f} {result['encirclements']:<15.2f} {result['Z']:<20.0f} {result['stability']:<10}")

这个实现虽然简化了，但包含了Nyquist判据的核心思想。在实际工程中，我们还需要考虑：

• 虚轴上的极点（需要做无穷小半圆绕行）
• 多圈包围的情况
• 频率从-∞到+∞的完整路径

3.3 处理有积分环节的系统

有积分环节的系统（即开环传递函数在原点有极点）在绘制Nyquist图时需要特殊处理，因为ω=0时幅值无穷大。control库的nyquist_plot()函数会自动处理这种情况，但了解原理很重要：

# 创建含有积分环节的系统
# G(s) = K/(s(s+1)(s+2))
K_values = [1, 5, 10]  # 不同增益

plt.figure(figsize=(15, 5))

for i, K in enumerate(K_values):
    num = [K]
    den = [1, 3, 2, 0]  # s(s+1)(s+2) = s^3 + 3s^2 + 2s
    sys = control.TransferFunction(num, den)
    
    plt.subplot(1, 3, i+1)
    
    # 绘制Nyquist图
    control.nyquist_plot(sys, omega_limits=(0.01, 100), 
                         arrow_width=0.015, arrow_color='blue')
    
    # 标记(-1, j0)点
    plt.plot(-1, 0, 'ro', markersize=10, markeredgewidth=2, 
             markerfacecolor='none', label='(-1, j0)')
    
    # 添加网格和坐标轴
    plt.grid(True, alpha=0.5)
    plt.axhline(y=0, color='k', linestyle='-', alpha=0.3, linewidth=0.5)
    plt.axvline(x=0, color='k', linestyle='-', alpha=0.3, linewidth=0.5)
    
    # 计算稳定裕度
    gm, pm, wg, wc = control.margin(sys)
    
    plt.title(f'K={K}\nPM={pm:.1f}°, GM={20*np.log10(gm):.1f}dB', fontsize=12)
    plt.xlabel('实部', fontsize=11)
    plt.ylabel('虚部', fontsize=11)
    plt.legend(loc='best', fontsize=9)
    plt.axis('equal')

plt.suptitle('含有积分环节系统的Nyquist图（不同增益K）', fontsize=14)
plt.tight_layout()
plt.show()

对于有积分环节的系统，Nyquist曲线在ω→0时会趋向无穷远。在实际绘图中，control库会从一个小正频率开始，避免数值问题。从图中我们可以清楚地看到：

• 增益K较小时，曲线不包围(-1, j0)点，系统稳定
• 增益K增大到一定程度，曲线开始包围(-1, j0)点，系统变得不稳定
• 通过Nyquist图可以直观地看出增益裕度

3.4 复杂系统的Nyquist分析

对于高阶系统或非最小相位系统，Nyquist图可能会比较复杂。这时我们需要更仔细地分析：

# 创建一个非最小相位系统示例
# G(s) = (1 - s)/((s+1)(s+2))
num = [-1, 1]  # 注意：1 - s = -s + 1，所以系数是[-1, 1]
den = [1, 3, 2]  # (s+1)(s+2) = s^2 + 3s + 2

sys_nonmin = control.TransferFunction(num, den)

# 绘制Nyquist图
plt.figure(figsize=(12, 10))

# 完整Nyquist图
plt.subplot(2, 2, 1)
control.nyquist_plot(sys_nonmin, omega_limits=(0.01, 100))
plt.plot(-1, 0, 'ro', markersize=10, markeredgewidth=2, 
         markerfacecolor='none', label='(-1, j0)')
plt.grid(True)
plt.title('非最小相位系统Nyquist图（完整）', fontsize=12)
plt.legend()

# 局部放大图（关注(-1, j0)点附近）
plt.subplot(2, 2, 2)
control.nyquist_plot(sys_nonmin, omega_limits=(0.01, 100))
plt.plot(-1, 0, 'ro', markersize=10, markeredgewidth=2, 
         markerfacecolor='none')
plt.grid(True)
plt.xlim([-3, 1])
plt.ylim([-2, 2])
plt.title('(-1, j0)点附近局部放大', fontsize=12)

# Bode图对比
plt.subplot(2, 2, 3)
control.bode(sys_nonmin, dB=True, deg=True, plot=True)
plt.suptitle('非最小相位系统Bode图', fontsize=12)

# 阶跃响应
plt.subplot(2, 2, 4)
t, y = control.step_response(control.feedback(sys_nonmin, 1))
plt.plot(t, y, 'b-', linewidth=2)
plt.grid(True)
plt.xlabel('时间 (s)', fontsize=11)
plt.ylabel('输出', fontsize=11)
plt.title('闭环系统阶跃响应', fontsize=12)

plt.tight_layout()
plt.show()

# 分析系统特性
print("=== 非最小相位系统分析 ===")
print(f"传递函数: {sys_nonmin}")
print(f"极点: {control.pole(sys_nonmin)}")
print(f"零点: {control.zero(sys_nonmin)}")

# 计算稳定裕度
gm, pm, wg, wc = control.margin(sys_nonmin)
print(f"\n稳定裕度分析:")
print(f"相位裕度: {pm:.2f} 度")
print(f"增益裕度: {20*np.log10(gm):.2f} dB")
print(f"相位穿越频率: {wg:.2f} rad/s")
print(f"截止频率: {wc:.2f} rad/s")

非最小相位系统（有右半平面零点）的Nyquist图通常会有一些反直觉的特征。通过多图对比分析，我们可以更全面地理解系统特性。

4. 实战应用：从理论到工程实践

掌握了Bode图和Nyquist图的绘制方法后，我们来看看如何将这些工具应用到实际工程问题中。我结合自己多年的项目经验，分享几个典型的应用场景。

4.1 控制系统设计与校正

假设我们需要设计一个位置伺服系统，要求相位裕度大于45度，增益裕度大于10dB。原始系统的开环传递函数为：

G(s) = 10/(s(s+1)(s+2))

让我们先分析原始系统的性能：

# 原始系统
num_orig = [10]
den_orig = [1, 3, 2, 0]  # s(s+1)(s+2)
sys_orig = control.TransferFunction(num_orig, den_orig)

# 计算原始系统性能
gm_orig, pm_orig, wg_orig, wc_orig = control.margin(sys_orig)
print("=== 原始系统性能 ===")
print(f"相位裕度: {pm_orig:.1f}° (要求>45°)")
print(f"增益裕度: {20*np.log10(gm_orig):.1f} dB (要求>10dB)")
print(f"截止频率: {wc_orig:.2f} rad/s")

# 设计超前校正器
# 超前校正器传递函数: Gc(s) = (1+αTs)/(1+Ts), α>1
# 我们需要增加约30度的相位裕度
pm_required = 45
pm_needed = pm_required - pm_orig + 10  # 增加10度裕量

# 计算超前校正器参数
phi_m = np.radians(pm_needed)  # 转换为弧度
alpha = (1 + np.sin(phi_m)) / (1 - np.sin(phi_m))
print(f"\n=== 超前校正器设计 ===")
print(f"需要增加的相位裕度: {pm_needed:.1f}°")
print(f"α值: {alpha:.3f}")

# 选择转折频率为原始系统的截止频率
T = 1 / (wc_orig * np.sqrt(alpha))
print(f"时间常数T: {T:.3f} s")
print(f"转折频率1/T: {1/T:.2f} rad/s")
print(f"转折频率1/(αT): {1/(alpha*T):.2f} rad/s")

# 创建校正器
num_c = [alpha*T, 1]
den_c = [T, 1]
sys_comp = control.TransferFunction(num_c, den_c)

# 校正后系统
sys_corrected = sys_comp * sys_orig

# 计算校正后性能
gm_corr, pm_corr, wg_corr, wc_corr = control.margin(sys_corrected)
print(f"\n=== 校正后系统性能 ===")
print(f"相位裕度: {pm_corr:.1f}°")
print(f"增益裕度: {20*np.log10(gm_corr):.1f} dB")
print(f"截止频率: {wc_corr:.2f} rad/s")

# 绘制对比图
plt.figure(figsize=(14, 10))

# Bode图对比
plt.subplot(2, 2, 1)
mag_orig, phase_orig, omega = control.bode(sys_orig, plot=False)
mag_corr, phase_corr, _ = control.bode(sys_corrected, omega, plot=False)

plt.semilogx(omega, 20*np.log10(mag_orig), 'b-', linewidth=2, label='原始系统')
plt.semilogx(omega, 20*np.log10(mag_corr), 'r--', linewidth=2, label='校正后系统')
plt.grid(True, which='both', alpha=0.5)
plt.ylabel('幅值 (dB)', fontsize=12)
plt.title('幅频特性对比', fontsize=13)
plt.legend()

plt.subplot(2, 2, 3)
plt.semilogx(omega, phase_orig, 'b-', linewidth=2, label='原始系统')
plt.semilogx(omega, phase_corr, 'r--', linewidth=2, label='校正后系统')
plt.grid(True, which='both', alpha=0.5)
plt.xlabel('频率 (rad/s)', fontsize=12)
plt.ylabel('相位 (度)', fontsize=12)
plt.title('相频特性对比', fontsize=13)
plt.legend()

# Nyquist图对比
plt.subplot(2, 2, 2)
control.nyquist_plot(sys_orig, omega_limits=(0.01, 100), 
                     label='原始系统', arrow_width=0.01)
control.nyquist_plot(sys_corrected, omega_limits=(0.01, 100), 
                     label='校正后系统', arrow_width=0.01)
plt.plot(-1, 0, 'ko', markersize=8, label='(-1, j0)')
plt.grid(True, alpha=0.5)
plt.axis('equal')
plt.title('Nyquist图对比', fontsize=13)
plt.legend()

# 阶跃响应对比
plt.subplot(2, 2, 4)
t_orig, y_orig = control.step_response(control.feedback(sys_orig, 1))
t_corr, y_corr = control.step_response(control.feedback(sys_corrected, 1))

plt.plot(t_orig, y_orig, 'b-', linewidth=2, label='原始系统')
plt.plot(t_corr, y_corr, 'r--', linewidth=2, label='校正后系统')
plt.grid(True, alpha=0.5)
plt.xlabel('时间 (s)', fontsize=12)
plt.ylabel('输出', fontsize=12)
plt.title('阶跃响应对比', fontsize=13)
plt.legend()

plt.suptitle('超前校正效果对比分析', fontsize=14)
plt.tight_layout()
plt.show()

这个例子展示了完整的控制系统设计流程。通过频域分析，我们可以：

1. 评估原始系统的性能缺陷
2. 设计合适的校正器
3. 验证校正效果
4. 通过时域响应确认性能改善

4.2 系统辨识与模型验证

在实际工程中，我们经常需要通过实验数据来辨识系统模型。频域方法在这方面特别有用：

# 模拟实验数据（在实际应用中，这些数据来自实际测量）
def generate_experimental_data(sys_true, noise_level=0.05):
    """生成带噪声的实验数据"""
    omega_exp = np.logspace(-1, 2, 50)  # 实验频率点
    
    # 计算真实频率响应
    mag_true, phase_true, _ = control.bode(sys_true, omega_exp, plot=False)
    
    # 添加测量噪声
    np.random.seed(42)  # 固定随机种子，确保可重复性
    mag_noise = mag_true * (1 + noise_level * np.random.randn(len(omega_exp)))
    phase_noise = phase_true + noise_level * 10 * np.random.randn(len(omega_exp))
    
    return omega_exp, mag_noise, phase_noise, mag_true, phase_true

# 真实系统（我们不知道，需要辨识）
sys_true = control.TransferFunction([5], [1, 2, 5, 1])  # 三阶系统

# 生成实验数据
omega_exp, mag_exp, phase_exp, mag_true, phase_true = generate_experimental_data(
    sys_true, noise_level=0.08
)

# 系统辨识：尝试用二阶模型拟合
def fit_transfer_function(omega, mag_db, phase_deg, order=2):
    """使用频域数据拟合传递函数"""
    from scipy import optimize
    
    # 将dB转换为线性值
    mag_linear = 10**(mag_db/20)
    
    # 构建频率响应数据
    freqs = 1j * omega
    response_data = mag_linear * np.exp(1j * np.radians(phase_deg))
    
    # 定义误差函数
    def error_func(params):
        """计算模型与实验数据的误差"""
        if order == 2:
            # 二阶系统: b0/(s^2 + a1*s + a0)
            b0 = params[0]
            a0, a1 = params[1], params[2]
            num = [b0]
            den = [1, a1, a0]
        elif order == 3:
            # 三阶系统: (b0)/(s^3 + a2*s^2 + a1*s + a0)
            b0 = params[0]
            a0, a1, a2 = params[1], params[2], params[3]
            num = [b0]
            den = [1, a2, a1, a0]
        else:
            raise ValueError("只支持2阶或3阶")
        
        sys = control.TransferFunction(num, den)
        response_model = sys(1j * omega)
        
        # 计算复数误差（幅值和相位）
        error = np.sum(np.abs(response_model - response_data)**2)
        return error
    
    # 初始猜测（根据经验）
    if order == 2:
        initial_guess = [5, 5, 2]  # [b0, a0, a1]
        bounds = [(0.1, 20), (0.1, 20), (0.1, 10)]
    else:
        initial_guess = [5, 1, 5, 2]  # [b0, a0, a1, a2]
        bounds = [(0.1, 20), (0.1, 10), (0.1, 20), (0.1, 10)]
    
    # 优化
    result = optimize.minimize(error_func, initial_guess, bounds=bounds, 
                               method='L-BFGS-B')
    
    # 提取最优参数
    if order == 2:
        b0, a0, a1 = result.x
        sys_fit = control.TransferFunction([b0], [1, a1, a0])
    else:
        b0, a0, a1, a2 = result.x
        sys_fit = control.TransferFunction([b0], [1, a2, a1, a0])
    
    return sys_fit, result.fun

# 尝试二阶和三阶拟合
sys_fit_2, error_2 = fit_transfer_function(omega_exp, 
                                           20*np.log10(mag_exp), 
                                           phase_exp, order=2)

sys_fit_3, error_3 = fit_transfer_function(omega_exp, 
                                           20*np.log10(mag_exp), 
                                           phase_exp, order=3)

print("=== 系统辨识结果 ===")
print(f"真实系统: {sys_true}")
print(f"\n二阶拟合系统: {sys_fit_2}")
print(f"拟合误差: {error_2:.4f}")
print(f"\n三阶拟合系统: {sys_fit_3}")
print(f"拟合误差: {error_3:.4f}")

# 绘制拟合结果对比
plt.figure(figsize=(14, 8))

# 幅频特性
plt.subplot(2, 2, 1)
# 实验数据
plt.semilogx(omega_exp, 20*np.log10(mag_exp), 'ko', 
             markersize=6, label='实验数据', alpha=0.7)
# 真实系统
mag_true_plot, phase_true_plot, _ = control.bode(sys_true, plot=False)
plt.semilogx(omega_exp, 20*np.log10(mag_true_plot), 'b-', 
             linewidth=2, label='真实系统')
# 拟合系统
mag_fit2, phase_fit2, _ = control.bode(sys_fit_2, omega_exp, plot=False)
mag_fit3, phase_fit3, _ = control.bode(sys_fit_3, omega_exp, plot=False)

plt.semilogx(omega_exp, 20*np.log10(mag_fit2), 'r--', 
             linewidth=2, label='二阶拟合')
plt.semilogx(omega_exp, 20*np.log10(mag_fit3), 'g-.', 
             linewidth=2, label='三阶拟合')

plt.grid(True, which='both', alpha=0.5)
plt.ylabel('幅值 (dB)', fontsize=12)
plt.title('幅频特性拟合对比', fontsize=13)
plt.legend()

# 相频特性
plt.subplot(2, 2, 3)
plt.semilogx(omega_exp, phase_exp, 'ko', 
             markersize=6, label='实验数据', alpha=0.7)
plt.semilogx(omega_exp, phase_true_plot, 'b-', 
             linewidth=2, label='真实系统')
plt.semilogx(omega_exp, phase_fit2, 'r--', 
             linewidth=2, label='二阶拟合')
plt.semilogx(omega_exp, phase_fit3, 'g-.', 
             linewidth=2, label='三阶拟合')

plt.grid(True, which='both', alpha=0.5)
plt.xlabel('频率 (rad/s)', fontsize=12)
plt.ylabel('相位 (度)', fontsize=12)
plt.title('相频特性拟合对比', fontsize=13)
plt.legend()

# Nyquist图对比
plt.subplot(2, 2, 2)
control.nyquist_plot(sys_true, omega_limits=(0.1, 100), 
                     label='真实系统', arrow_width=0.01)
control.nyquist_plot(sys_fit_2, omega_limits=(0.1, 100), 
                     label='二阶拟合', arrow_width=0.01)
control.nyquist_plot(sys_fit_3, omega_limits=(0.1, 100), 
                     label='三阶拟合', arrow_width=0.01)
plt.plot(-1, 0, 'ko', markersize=8)
plt.grid(True, alpha=0.5)
plt.axis('equal')
plt.title('Nyquist图对比', fontsize=13)
plt.legend()

# 阶跃响应对比
plt.subplot(2, 2, 4)
t = np.linspace(0, 10, 1000)
_, y_true = control.step_response(sys_true, t)
_, y_fit2 = control.step_response(sys_fit_2, t)
_, y_fit3 = control.step_response(sys_fit_3, t)

plt.plot(t, y_true, 'b-', linewidth=2, label='真实系统')
plt.plot(t, y_fit2, 'r--', linewidth=2, label='二阶拟合')
plt.plot(t, y_fit3, 'g-.', linewidth=2, label='三阶拟合')
plt.grid(True, alpha=0.5)
plt.xlabel('时间 (s)', fontsize=12)
plt.ylabel('输出', fontsize=12)
plt.title('阶跃响应对比', fontsize=13)
plt.legend()

plt.suptitle('基于频域数据的系统辨识', fontsize=14)
plt.tight_layout()
plt.show()

这个例子展示了如何利用频域数据进行系统辨识。在实际工程中，我们通常：

1. 通过扫频实验获取频率响应数据
2. 选择合适的模型结构（阶次）
3. 用优化算法拟合模型参数
4. 验证模型在时域和频域的准确性

4.3 鲁棒性分析与灵敏度函数

在现代控制系统中，鲁棒性是一个重要考量。我们可以通过频域分析来评估系统对参数变化的敏感度：

# 研究参数变化对系统性能的影响
def analyze_parameter_variation(nominal_sys, param_name, param_values):
    """分析参数变化对系统性能的影响"""
    
    results = []
    
    for value in param_values:
        # 根据参数名修改系统
        if param_name == '增益K':
            sys = value * nominal_sys
        elif param_name == '时间常数T':
            # 假设系统为 G(s) = K/(Ts+1)
            K = 1
            sys = control.TransferFunction([K], [value, 1])
        elif param_name == '阻尼比ζ':
            # 假设系统为 G(s) = ω_n^2/(s^2 + 2ζω_n s + ω_n^2)
            omega_n = 1
            sys = control.TransferFunction([omega_n**2], 
                                          [1, 2*value*omega_n, omega_n**2])
        else:
            continue
        
        # 计算性能指标
        gm, pm, wg, wc = control.margin(sys)
        
        # 计算时域指标
        t, y = control.step_response(control.feedback(sys, 1))
        
        # 超调量
        y_ss = y[-1]
        overshoot = (np.max(y) - y_ss) / y_ss * 100 if y_ss != 0 else 0
        
        # 调节时间（2%误差带）
        settling_idx = np.where(np.abs(y - y_ss) <= 0.02 * np.abs(y_ss))[0]
        settling_time = t[settling_idx[0]] if len(settling_idx) > 0 else t[-1]
        
        # 上升时间（10%到90%）
        y_min, y_max = np.min(y), np.max(y)
        y_10 = y_min + 0.1 * (y_max - y_min)
        y_90 = y_min + 0.9 * (y_max - y_min)
        
        idx_10 = np.where(y >= y_10)[0][0]
        idx_90 = np.where(y >= y_90)[0][0]
        rise_time = t[idx_90] - t[idx_10]
        
        results.append({
            '参数值': value,
            '相位裕度': pm,
            '增益裕度_db': 20*np.log10(gm),
            '截止频率': wc,
            '超调量': overshoot,
            '调节时间': settling_time,
            '上升时间': rise_time
        })
    
    return results

# 示例：分析增益变化对系统性能的影响
nominal_sys = control.TransferFunction([10], [1, 3, 2, 0])
K_values = np.linspace(1, 20, 10)  # 增益从1到20

results = analyze_parameter_variation(nominal_sys, '增益K', K_values)

# 创建性能指标表格
print("=== 增益变化对系统性能的影响 ===")
print("-" * 90)
print(f"{'增益K':<8} {'相位裕度(°)':<12} {'增益裕度(dB)':<14} {'截止频率':<10} {'超调量(%)':<10} {'调节时间(s)':<12} {'上升时间(s)':<10}")
print("-" * 90)

for r in results:
    print(f"{r['参数值']:<8.2f} {r['相位裕度']:<12.2f} {r['增益裕度_db']:<14.2f} "
          f"{r['截止频率']:<10.2f} {r['超调量']:<10.2f} {r['调节时间']:<12.2f} {r['上升时间']:<10.2f}")

# 可视化分析
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

# 提取数据
K_vals = [r['参数值'] for r in results]
PM_vals = [r['相位裕度'] for r in results]
GM_vals = [r['增益裕度_db'] for r in results]
wc_vals = [r['截止频率'] for r in results]
OS_vals = [r['超调量'] for r in results]
Ts_vals = [r['调节时间'] for r in results]
Tr_vals = [r['上升时间'] for r in results]

# 相位裕度 vs 增益
axes[0, 0].plot(K_vals, PM_vals, 'bo-', linewidth=2, markersize=8)
axes[0, 0].axhline(y=45, color='r', linestyle='--', alpha=0.7, label='目标PM=45°')
axes[0, 0].grid(True, alpha=0.5)
axes[0, 0].set_xlabel('增益K', fontsize=12)
axes[0, 0].set_ylabel('相位裕度 (°)', fontsize=12)
axes[0, 0].set_title('相位裕度 vs 增益', fontsize=13)
axes[0, 0].legend()

# 增益裕度 vs 增益
axes[0, 1].plot(K_vals, GM_vals, 'go-', linewidth=2, markersize=8)
axes[0, 1].axhline(y=10, color='r', linestyle='--', alpha=0.7, label='目标GM=10dB')
axes[0, 1].grid(True, alpha=0.5)
axes[0, 1].set_xlabel('增益K', fontsize=12)
axes[0, 1].set_ylabel('增益裕度 (dB)', fontsize=12)
axes[0, 1].set_title('增益裕度 vs 增益', fontsize=13)
axes[0, 1].legend()

# 截止频率 vs 增益
axes[0, 2].plot(K_vals, wc_vals, 'ro-', linewidth=2, markersize=8)
axes[0, 2].grid(True, alpha=0.5)
axes[0, 2].set_xlabel('增益K', fontsize=12)
axes[0, 2].set_ylabel('截止频率 (rad/s)', fontsize=12)
axes[0, 2].set_title('截止频率 vs 增益', fontsize=13)

# 超调量 vs 增益
axes[1, 0].plot(K_vals, OS_vals, 'mo-', linewidth=2, markersize=8)
axes[1, 0].axhline(y=20, color='r', linestyle='--', alpha=0.7, label='目标OS<20%')
axes[1, 0].grid(True, alpha=0.5)
axes[1, 0].set_xlabel('增益K', fontsize=12)
axes[1, 0].set_ylabel('超调量 (%)', fontsize=12)
axes[1, 0].set_title('超调量 vs 增益', fontsize=13)
axes[1, 0].legend()

# 调节时间 vs 增益
axes[1, 1].plot(K_vals, Ts_vals, 'co-', linewidth=2, markersize=8)
axes[1, 1].axhline(y=2, color='r', linestyle='--', alpha=0.7, label='目标Ts<2s')
axes[1, 1].grid(True, alpha=0.5)
axes[1, 1].set_xlabel('增益K', fontsize=12)
axes[1, 1].set_ylabel('调节时间 (s)', fontsize=12)
axes[1, 1].set_title('调节时间 vs 增益', fontsize=13)
axes[1, 1].legend()

# 上升时间 vs 增益
axes[1, 2].plot(K_vals, Tr_vals, 'yo-', linewidth=2, markersize=8)
axes[1, 2].grid(True, alpha=0.5)
axes[1, 2].set_xlabel('增益K', fontsize=12)
axes[1, 2].set_ylabel('上升时间 (s)', fontsize=12)
axes[1, 2].set_title('上升时间 vs 增益', fontsize=13)

plt.suptitle('系统性能对增益变化的敏感性分析', fontsize=14)
plt.tight_layout()
plt.show()

# 找到满足所有性能要求的增益范围
valid_K = []
for r in results:
    if (r['相位裕度'] >= 45 and 
        r['增益裕度_db'] >= 10 and 
        r['超调量'] <= 20 and 
        r['调节时间'] <= 2):
        valid_K.append(r['参数值'])

if valid_K:
    print(f"\n满足所有性能要求的增益范围: {min(valid_K):.2f} ≤ K ≤ {max(valid_K):.2f}")
else:
    print("\n没有增益值能同时满足所有性能要求")
    print("可能需要更复杂的校正策略")

这种敏感性分析在实际工程中非常有用。它帮助我们：

1. 理解参数变化对系统性能的影响
2. 确定参数的允许变化范围
3. 评估系统的鲁棒性
4. 为控制器设计提供指导

4.4 批量处理与自动化报告生成

在实际工程项目中，我们经常需要分析多个系统或同一系统的多个变体。这时，自动化处理就变得非常重要：

import pandas as pd
from datetime import datetime

def analyze_multiple_systems(systems_dict):
    """批量分析多个系统"""
    
    results = []
    
    for name, sys in systems_dict.items():
        # 计算频域指标
        try:
            gm, pm, wg, wc = control.margin(sys)
            gm_db = 20 * np.log10(gm)
        except:
            gm, pm, wg, wc = np.nan, np.nan, np.nan, np.nan
            gm_db = np.nan
        
        # 计算时域指标
        try:
            t, y = control.step_response(control.feedback(sys, 1), T=np.linspace(0, 10, 1000))
            
            # 稳态值
            y_ss = y[-1]
            
            # 超调量
            overshoot = (np.max(y) - y_ss) / y_ss * 100 if y_ss != 0 else 0
            
            # 调节时间（2%误差带）
            settling_idx = np.where(np.abs(y - y_ss) <= 0.02 * np.abs(y_ss))[0]
            settling_time = t[settling_idx[0]] if len(settling_idx) > 0 else t[-1]
            
            # 上升时间（10%到90%）
            y_min, y_max = np.min(y), np.max(y)
            y_10 = y_min + 0.1 * (y_max - y_min)
            y_90 = y_min + 0.9 * (y_max - y_min)
            
            idx_10 = np.where(y >= y_10)[0][0] if len(np.where(y >= y_10)[0]) > 0 else 0
            idx_90 = np.where(y >= y_90)[0][0] if len(np.where(y >= y_90)[0]) > 0 else 0
            rise_time = t[idx_90] - t[idx_10] if idx_90 > idx_10 else np.nan
            
        except:
            overshoot = np.nan
            settling_time = np.nan
            rise_time = np.nan
        
        # 计算带宽
        try:
            # 计算闭环频率响应
            sys_closed = control.feedback(sys, 1)
            omega = np.logspace(-2, 2, 1000)
            mag, phase, omega = control.bode(sys_closed, omega, plot=False)
            
            # 零频率增益
            mag_0 = np.abs(control.evalfr(sys_closed, 0))
            mag_0_db = 20 * np.log10(mag_0)
            
            # 带宽频率（-3dB点）
            bandwidth_idx = np.where(20*np.log10(mag) <= mag_0_db - 3)[0]
            bandwidth = omega[bandwidth_idx[0]] if len(bandwidth_idx) > 0 else np.nan
            
        except:
            bandwidth = np.nan
        
        # 存储结果
        results.append({
            '系统名称': name,
            '传递函数': str(sys),
            '相位裕度(°)': f"{pm:.2f}" if not np.isnan(pm) else "N/A",
            '增益裕度(dB)': f"{gm_db:.2f}" if not np.isnan(gm_db) else "N/A",
            '截止频率(rad/s)': f"{wc:.2f}" if not np.isnan(wc) else "N/A",
            '带宽(rad/s)': f"{bandwidth:.2f}" if not np.isnan(bandwidth) else "N/A",
            '超调量(%)': f"{overshoot:.2f}" if not np.isnan(overshoot) else "N/A",
            '调节时间(s)': f"{settling_time:.2f}" if not np.isnan(settling_time) else "N/A",
            '上升时间(s)': f"{rise_time:.2f}" if not np.isnan(rise_time) else "N/A",
            '稳定性': '稳定' if pm > 0 and gm_db > 0 else '不稳定'
        })
    
    return pd.DataFrame(results)

# 定义多个测试系统
test_systems = {
    '二阶欠阻尼': control.TransferFunction([1], [1, 0.8, 1]),
    '二阶临界阻尼': control.TransferFunction([1], [1, 2, 1]),
    '二阶过阻尼': control.TransferFunction([1], [1, 3, 1]),
    '积分环节': control.TransferFunction([1], [1, 0]),
    '一阶惯性': control.TransferFunction([1], [2, 1]),
    '不稳定系统': control.TransferFunction([1], [1, -0.5, 1]),
    '三阶系统': control.TransferFunction([10], [1, 3, 2, 0]),
    '带零点系统': control.TransferFunction([1, 1], [1, 2, 1]),
}

# 执行批量分析
df_results = analyze_multiple_systems(test_systems)

# 显示结果表格
print("=== 多系统性能分析报告 ===")
print(f"生成时间: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')}")
print("=" * 120)
print(df_results.to_string(index=False))
print("=" * 120)

# 保存到CSV文件
df_results.to_csv('system_analysis_report.csv', index=False, encoding='utf-8-sig')
print("\n分析报告已保存到: system_analysis_report.csv")

# 可视化比较
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 选择几个典型系统进行详细比较
selected_systems = ['二阶欠阻尼', '二阶临界阻尼', '二阶过阻尼', '三阶系统']
colors = ['b', 'g', 'r', 'm']

# Bode图比较
ax1 = axes[0, 0]
for i, name in enumerate(selected_systems):
    sys = test_systems[name]
    mag, phase, omega = control.bode(sys, plot=False)
    ax1.semilogx(omega, 20*np.log10(mag), color=colors[i], 
                 linewidth=2, label=name)
ax1.grid(True, which='both', alpha=0.5)
ax1.set_ylabel('幅值 (dB)', fontsize=12)
ax1.set_title('Bode图幅频特性比较', fontsize=13)
ax1.legend()

ax2 = axes[1, 0]
for i, name in enumerate(selected_systems):
    sys = test_systems[name]
    mag, phase, omega = control.bode(sys, plot=False)
    ax2.semilogx(omega, phase, color=colors[i], 
                 linewidth=2, label=name)
ax2.grid(True, which='both', alpha=0.5)
ax2.set_xlabel('频率 (rad/s)', fontsize=12)
ax2.set_ylabel('相位 (度)', fontsize=12)
ax2.set_title('Bode图相频特性比较', fontsize=13)
ax2.legend()

# Nyquist图比较
ax3 = axes[0, 1]
for i, name in enumerate(selected_systems):
    sys = test_systems[name]
    # 计算频率响应
    omega = np.logspace(-2, 2, 1000)
    response = sys(1j * omega)
    ax3.plot(np.real(response), np.imag(response), 
             color=colors[i], linewidth=2, label=name)
ax3.plot(-1, 0, 'ko', markersize=8, label='(-1, j0)')
ax3.grid(True, alpha=0.5)
ax3.axhline(y=0, color='k', linestyle='-', alpha=0.3, linewidth=0.5)
ax3.axvline(x=0, color='k', linestyle='-', alpha=0.3, linewidth=0.5)
ax3.set_xlabel('实部', fontsize=12)
ax3.set_ylabel('虚部', fontsize=12)
ax3.set_title('Nyquist图比较', fontsize=13)
ax3.legend()
ax3.axis('equal')

# 阶跃响应比较
ax4 = axes[1, 1]
t = np.linspace(0, 10, 1000)
for i, name in enumerate(selected_systems):
    sys = test_systems[name]
    try:
        t, y = control.step_response(control.feedback(sys, 1), t)
        ax4.plot(t, y, color=colors[i], linewidth=2, label=name)
    except:
        continue
ax4.grid(True, alpha=0.5)
ax4.set_xlabel('时间 (s)', fontsize=12)
ax4.set_ylabel('输出', fontsize=12)
ax4.set_title('阶跃响应比较', fontsize=13)
ax4.legend()

plt.suptitle('多系统性能对比分析', fontsize=14)
plt.tight_layout()
plt.show()

# 生成性能雷达图
def create_radar_chart(system_names, metrics_df):
    """创建性能雷达图"""
    
    # 选择要展示的指标
    metrics = ['相位裕度(°)', '增益裕度(dB)', '带宽(rad/s)', '超调量(%)', '调节时间(s)']
    
    # 准备数据
    categories = metrics
    N = len(categories)
    
    # 创建子图
    fig, axes = plt.subplots(2, 2, figsize=(12, 10), 
                            subplot_kw=dict(projection='polar'))
    axes = axes.flatten()
    
    for idx, sys_name in enumerate(system_names[:4]):  # 只显示前4个系统
        ax = axes[idx]
        
        # 获取数据
        row = metrics_df[metrics_df['系统名称'] == sys_name].iloc[0]
        
        values = []
        for metric in metrics:
            val_str = row[metric]
            if val_str != 'N/A':
                values.append(float(val_str.split()[0]))  # 提取数值部分
            else:
                values.append(0)
        
        # 归一化处理（使所有指标在0-1范围内）
        # 这里需要根据实际情况调整归一化方式
        normalized_values = []
        for i, val in enumerate(values):
            if metrics[i] in ['超调量(%)', '调节时间(s)']:
                # 这些指标越小越好，所以用倒数
                normalized = 1.0 / (1 + val/100) if val > 0 else 0
            else:
                # 这些指标越大越好
                normalized = val / 100 if val <= 100 else 1.0
            normalized_values.append(normalized)
        
        # 完成循环
        values_radar = normalized_values + normalized_values[:1]
        
        # 计算角度
        angles = [n / float(N) * 2 * np.pi for n in range(N)]
        angles += angles[:1]
        
        # 绘制雷达图
        ax.plot(angles, values_radar, 'o-', linewidth=2, label=sys_name)
        ax.fill(angles, values_radar, alpha=0.25)
        
        # 设置标签
        ax.set_xticks(angles[:-1])
        ax.set_xticklabels(categories, fontsize=9)
        ax.set_ylim(0, 1)
        ax.set_title(sys_name, fontsize=11, pad=20)
    
    plt.suptitle('系统性能雷达图（归一化比较）', fontsize=14)
    plt.tight_layout()
    plt.show()

# 创建雷达图
create_radar_chart(selected_systems, df_results)

这种批量分析方法特别适合：

• 比较多个设计方案
• 参数优化研究
• 系统选型评估
• 生成项目报告

通过自动化分析，我们可以快速评估大量系统变体，找到最优设计方案，并生成专业的分析报告。