企业官网建设流程全解析

从频响曲线到奈奎斯特图：用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 不同零点配置的影响分析

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

1. T3 > T2 > T1（零点转折频率最小）
2. T3 < T2 < T1（零点转折频率最大）
3. 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. 实用技巧与常见问题解决

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

1. 高频区域分辨率不足：

   • 增加高频区域的采样点

   omega = np.logspace(-1, 3, 2000) # 扩展到更高频率

2. 曲线不闭合：

   • 确保包含足够低的频率点
   • 手动添加s=0和s=∞的点

3. 数值不稳定：

   • 使用更高精度的计算

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

4. 多极点系统处理：

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

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

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

在实际项目中，我发现将奈奎斯特图与伯德图结合分析特别有效。先用伯德图快速了解系统大致特性，再用奈奎斯特图精确判断稳定性边界。对于复杂的多环路系统，可以分段绘制奈奎斯特图来简化分析。