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

从频响曲线到奈奎斯特图用Python实现控制系统稳定性可视化分析在控制系统的设计与分析中奈奎斯特图Nyquist Plot是一种强大的工具它能直观展示系统在不同频率下的响应特性。与伯德图Bode Plot相比奈奎斯特图将幅频和相频特性整合在复平面上使得系统稳定性的判断更加直观。本文将带你用PythonNumPyMatplotlib从零开始构建奈奎斯特图并通过实际代码演示如何分析不同系统配置下的稳定性特征。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, labelPositive Frequencies) plt.plot(H.real, -H.imag, r--, labelNegative 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]), textcoordsoffset points, xytext(10,5), hacenter) 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零点转折频率居中我们可以用以下代码生成对比图# 定义三种情况的传递函数 # 情况1T3 T2 T1 num1 [T3, 1] den1 [T1*T2, T1T2, 1, 0] sys1 signal.TransferFunction(num1, den1) # 情况2T3 T2 T1 num2 [T3, 1] den2 [T1*T2, T1T2, 1, 0] sys2 signal.TransferFunction(num2, den2) # 情况3T2 T3 T1 num3 [T3, 1] den3 [T1*T2, T1T2, 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, labelCase 1: T3 T2 T1) plt.plot(H2.real, H2.imag, r, labelCase 2: T3 T2 T1) plt.plot(H3.real, H3.imag, g, labelCase 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(fSystem is stable: {is_stable(sys, omega)})该函数首先计算频率响应然后通过分析相位变化来计算绕数最后根据奈奎斯特判据给出稳定性结论。5. 实用技巧与常见问题解决在实际应用中可能会遇到以下问题及解决方案高频区域分辨率不足增加高频区域的采样点omega np.logspace(-1, 3, 2000) # 扩展到更高频率曲线不闭合确保包含足够低的频率点手动添加s0和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第三象限混合特性中等稳定在实际项目中我发现将奈奎斯特图与伯德图结合分析特别有效。先用伯德图快速了解系统大致特性再用奈奎斯特图精确判断稳定性边界。对于复杂的多环路系统可以分段绘制奈奎斯特图来简化分析。