从频响曲线到奈奎斯特图：手把手教你用Python（NumPy+Matplotlib）绘制与分析控制系统稳定性

Jonna轩姐

253人浏览 · 2026-06-06 11:12:30

Jonna轩姐 · 2026-06-06 11:12:30 发布

从频响曲线到奈奎斯特图：用Python实现控制系统稳定性可视化分析

在控制系统的设计与分析中，奈奎斯特图（Nyquist Plot）是一种强大的工具，它能直观展示系统在不同频率下的响应特性。与伯德图（Bode Plot）相比，奈奎斯特图将幅频和相频特性整合在复平面上，使得系统稳定性的判断更加直观。本文将带你用Python（NumPy+Matplotlib）从零开始构建奈奎斯特图，并通过实际代码演示如何分析不同系统配置下的稳定性特征。

1. 理论基础与准备工作

奈奎斯特稳定性判据的核心思想是：通过分析开环传递函数G(jω)H(jω)在复平面上的轨迹，判断闭环系统的稳定性。具体来说，我们需要观察奈奎斯特曲线是否包围(-1, j0)点。

在开始编程前，需要准备以下Python库：

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

对于传递函数，我们通常表示为：

$$ G(s) = \frac{N(s)}{D(s)} = \frac{b_ms^m + b_{m-1}s^{m-1} + ... + b_0}{a_ns^n + a_{n-1}s^{n-1} + ... + a_0} $$

在Python中，我们可以使用 scipy.signal 模块的 TransferFunction 类来表示：

# 示例：创建传递函数 G(s) = 1/(s^2 + 2s + 1)
num = [1]  # 分子系数
den = [1, 2, 1]  # 分母系数
sys = signal.TransferFunction(num, den)

2. 频率响应计算与基本奈奎斯特图绘制

要绘制奈奎斯特图，首先需要计算系统在不同频率下的响应。NumPy的 logspace 函数非常适合生成频率点：

# 生成对数均匀分布的频率点（0.1到100 rad/s之间）
omega = np.logspace(-1, 2, 1000)

# 计算频率响应
w, H = signal.freqresp(sys, omega)

得到的 H 是一个复数数组，表示系统在各个频率点上的响应。我们可以直接绘制其实部和虚部：

plt.figure(figsize=(8, 6))
plt.plot(H.real, H.imag, 'b')
plt.plot(H.real, -H.imag, 'r--')  # 负频率部分
plt.grid(True)
plt.xlabel('Real')
plt.ylabel('Imaginary')
plt.title('Basic Nyquist Plot')
plt.show()

注意：奈奎斯特图应当包含正负频率两部分，因此我们同时绘制了H和其共轭。

3. 高级奈奎斯特图特性实现

3.1 自动标记关键频率点

在实际工程分析中，我们常需要关注转折频率等关键点。以下代码实现了自动标记功能：

def plot_nyquist_with_critical_points(sys, omega):
    w, H = signal.freqresp(sys, omega)
    
    plt.figure(figsize=(10, 8))
    plt.plot(H.real, H.imag, 'b', label='Positive Frequencies')
    plt.plot(H.real, -H.imag, 'r--', label='Negative Frequencies')
    
    # 标记关键频率点
    critical_points = [0.1, 1, 10]  # 示例关键频率
    for wp in critical_points:
        idx = np.argmin(np.abs(w - wp))
        plt.plot(H.real[idx], H.imag[idx], 'ko')
        plt.annotate(f'ω={wp:.1f}', 
                    (H.real[idx], H.imag[idx]),
                    textcoords="offset points",
                    xytext=(10,5),
                    ha='center')
    
    plt.grid(True)
    plt.legend()
    plt.xlabel('Real')
    plt.ylabel('Imaginary')
    plt.title('Nyquist Plot with Critical Points')
    plt.show()

# 使用示例
plot_nyquist_with_critical_points(sys, omega)

3.2 不同零点配置的影响分析

考虑以下三种情况的传递函数：

T3 > T2 > T1 （零点转折频率最小）
T3 < T2 < T1 （零点转折频率最大）
T2 > T3 > T1 （零点转折频率居中）

我们可以用以下代码生成对比图：

# 定义三种情况的传递函数
# 情况1：T3 > T2 > T1
num1 = [T3, 1]
den1 = [T1*T2, T1+T2, 1, 0]
sys1 = signal.TransferFunction(num1, den1)

# 情况2：T3 < T2 < T1
num2 = [T3, 1]
den2 = [T1*T2, T1+T2, 1, 0]
sys2 = signal.TransferFunction(num2, den2)

# 情况3：T2 > T3 > T1
num3 = [T3, 1]
den3 = [T1*T2, T1+T2, 1, 0]
sys3 = signal.TransferFunction(num3, den3)

# 计算频率响应
w1, H1 = signal.freqresp(sys1, omega)
w2, H2 = signal.freqresp(sys2, omega)
w3, H3 = signal.freqresp(sys3, omega)

# 绘制对比图
plt.figure(figsize=(12, 10))
plt.plot(H1.real, H1.imag, 'b', label='Case 1: T3 > T2 > T1')
plt.plot(H2.real, H2.imag, 'r', label='Case 2: T3 < T2 < T1')
plt.plot(H3.real, H3.imag, 'g', label='Case 3: T2 > T3 > T1')
plt.grid(True)
plt.legend()
plt.xlabel('Real')
plt.ylabel('Imaginary')
plt.title('Nyquist Plot Comparison with Different Zero Placements')
plt.show()

4. 稳定性自动判断算法实现

奈奎斯特稳定性判据的关键是判断曲线是否包围(-1, j0)点。我们可以实现一个简单的算法来自动完成这一判断：

def is_stable(sys, omega):
    w, H = signal.freqresp(sys, omega)
    H_full = np.concatenate([H, H.conj()[-2:0:-1]])  # 闭合曲线
    
    # 计算绕数
    angle_changes = np.diff(np.angle(H_full + 1 + 0j))
    angle_changes = (angle_changes + np.pi) % (2*np.pi) - np.pi
    winding_number = np.sum(angle_changes) / (2*np.pi)
    
    # 根据奈奎斯特判据判断稳定性
    poles = sys.poles
    unstable_open_loop_poles = sum(p.real > 0 for p in poles)
    return winding_number == -unstable_open_loop_poles / 2

# 使用示例
print(f"System is stable: {is_stable(sys, omega)}")

该函数首先计算频率响应，然后通过分析相位变化来计算绕数，最后根据奈奎斯特判据给出稳定性结论。

5. 实用技巧与常见问题解决

在实际应用中，可能会遇到以下问题及解决方案：

高频区域分辨率不足 ：
- 增加高频区域的采样点
```
omega = np.logspace(-1, 3, 2000)  # 扩展到更高频率
```
曲线不闭合 ：
- 确保包含足够低的频率点
- 手动添加s=0和s=∞的点

数值不稳定 ：

使用更高精度的计算

omega = np.logspace(-1, 2, 1000).astype(np.float64)

多极点系统处理 ：

对于在原点有多个极点的系统，需要采用小半圆绕行

# 创建绕行路径
epsilon = 1e-5
theta = np.linspace(-np.pi/2, np.pi/2, 100)
s_small = epsilon * np.exp(1j * theta)

以下表格总结了不同系统配置下的奈奎斯特图特征：

系统配置	起点象限	高频渐近线	稳定性特征
无零点	第三象限	沿虚轴接近	取决于增益
T3 > T2 > T1	第四象限	沿虚轴接近	更易稳定
T3 < T2 < T1	第三象限	沿实轴接近	可能振荡
T2 > T3 > T1	第三象限	混合特性	中等稳定