奈奎斯特曲线Nyquist plot是我做控制系统分析和现场调试时用得最多的工具之一。它比Bode图更直观一眼就能看出系统在临界点(-1, j0)附近的表现但也是让很多初学者头疼的东西频响、复数轨迹、稳定性判据混在一起稍不留神就画错方向、搞混坐标。这篇博文我想从实战角度出发把绘制奈奎斯特曲线的三种方法完整梳理一遍手工解析法、数值计算法、由Bode图转换法。三种方法各有适用场景不存在谁完全替代谁你把它们都吃透后面看任何教材里的奈奎斯特图都会觉得清清楚楚。这篇文章适合正在学自动控制原理的学生也适合刚入职做控制系统调试验证、需要快速判断稳定性的工程师。1. 画之前先把奈奎斯特图的基本逻辑搞透1.1 它到底在画什么奈奎斯特图的本质很简单把系统的开环传递函数 (G(s)) 中的 (s) 换成 (j\omega)得到随频率变化的复数 (G(j\omega))然后以实部为横轴、虚部为纵轴把 (\omega) 从 (0) 扫到 (\infty) 对应的所有复数点连成一条轨迹。很多人第一次看到这个概念会懵明明叫“频率响应”Bode图分成幅频和相频两条曲线奈奎斯特图却把它们合并成了一个点在复平面上移动。你可以这样理解Bode图把幅值和相角分开讲而奈奎斯特图把两者打包成一个极坐标点。一个点包含两项信息所以图上的横纵坐标不是频率而是复数的实部和虚部。举个例子最简单的惯性环节[ G(s) \frac{K}{Ts1} ]代入 (s j\omega)[ G(j\omega) \frac{K}{1j\omega T} ]把分母实数化[ G(j\omega) \frac{K(1-j\omega T)}{1\omega^2 T^2} ]实部是 (\frac{K}{1\omega^2 T^2})虚部是 (-\frac{K\omega T}{1\omega^2 T^2})。当 (\omega) 从 (0) 变到 (\infty) 时这条轨迹是一个圆心在 ((K/2, 0))、半径为 (K/2) 的下半圆。低频对应最右侧点高频逐渐朝原点靠拢。这就是奈奎斯特图的经典形态之一。1.2 为什么大家总是盯着(-1, j0)这个点稳定性是自动控制绕不开的话题。闭环系统的开环传函 (G(s)H(s)) 的奈奎斯特图和闭环特征方程之间有一条极其重要的数学联系也就是奈奎斯特判据。判据的常见表述是[ Z P N ]其中 (P) 是开环传递函数在右半 (s) 平面的极点数(N) 是当 (\omega) 从 (-\infty) 扫到 (\infty) 时奈奎斯特曲线顺时针绕过 ((-1, j0)) 点的净圈数(Z) 是闭环系统在右半平面的极点数。闭环稳定要求 (Z 0)。这个判据刚接触时容易写成死记硬背其实背后的逻辑并不难理解反馈系统闭环后稳定性发生变化的根源是特征方程根的分布而奈奎斯特曲线绕 ((-1, j0)) 的圈数恰好把开环频率特性和闭环特征根在右半平面的个数关联起来。你不需要每画一次图就推一遍围线积分但一定要记得这个结论是整个奈奎斯特方法的命根子。实际工程中大部分开环系统本身是稳定的也就是 (P0)这时候问题就简化成看曲线是否包围 ((-1, j0)) 点。不包围稳定包围不稳定。判断包围这件事比分析一堆代数特征方程要直观得多。1.3 画图前必须确认的三件事我见过太多人在画图时翻了车原因不是不会算而是没确认这三件事就急着上工具。第一开环传递函数是否已经是标准形式。建议先用零极点或时间常数的形式把它写清楚分母上有没有积分环节、分子上是超前还是滞后直接决定曲线的起点和终点方向。第二是否存在 (s0) 处的开环极点。如果系统类型是包含积分环节的例如 (\frac{K}{s(Ts1)})那么 (\omega \to 0) 时幅值趋向无穷大曲线起点在无穷远处。为了严格应用奈奎斯特判据需要补画一条从正频率零起点绕到负频率零起点的“无穷大半圆”否则包围圈数的判断会出错。第三确定需要的频率范围。稳定性主要看幅值接近 1 和相角接近 -180° 附近的中频段但手工草图通常要从低频一路画到高频才能看到全局走势数值计算则要单独选频率点选得太窄会漏掉关键的穿越点。2. 方法一手工解析法——先把原理吃透再谈工具2.1 核心流程代入sjω分离实频特性与虚频特性手工绘制的第一步永远是写出 (G(j\omega)) 的实部和虚部表达式。这一步熟练以后可以省略但在学习阶段绝对不能跳。标准套路有三步第一步写出开环传递函数 (G(s))然后把 (sj\omega) 代入。 第二步把分母中的虚数消掉。方法就是乘上分母的共轭把分母变成实数的模平方。 第三步整理成 (G(j\omega)R(\omega)jI(\omega)) 的形式得到实部 (R(\omega)) 和虚部 (I(\omega))。以 (G(s)\frac{1}{s(s1)}) 为例代入 (sj\omega)[ G(j\omega) \frac{1}{j\omega(1j\omega)} ]分母展开为 (j\omega - \omega^2)乘上共轭 ((-j\omega - \omega^2))整理后[ G(j\omega) \frac{-1}{1\omega^2} - j\frac{1}{\omega(1\omega^2)} ]于是[ R(\omega) -\frac{1}{1\omega^2},\quad I(\omega) -\frac{1}{\omega(1\omega^2)} ]注意一点这里更容易出错的是实部和虚部的符号。尤其在分母含有 (j\omega) 项时整理过程别跳步。算完后可以用一个具体的频率点验证比如 (\omega1) 时原式直接算一遍复数结果看跟你整理出的实部虚部是否一致。我平时哪怕用软件算也习惯保留这种手动校验的习惯。2.2 特殊点、渐近线、穿越实轴点手绘最关键的三步有了实部虚部表达式接下来不需要真的描几百个点抓住几个特殊点就能把曲线形状定个八九不离十。第一步是低频起点。低频指 (\omega \to 0^)。对 (\frac{1}{s(s1)})从上面的式子看实部趋向 (-1)虚部趋向负无穷大。也就是说曲线从实部约 (-1) 的地方出发但纵坐标在很远的负方向。严格讲是起点位于第三象限的无穷远处。第二步是高频终点。(\omega \to \infty) 时实部趋向 (0^-)虚部也趋向 (0^-)。分子分母阶次差是 2最终相角趋向 (-180^\circ)所以曲线从第三象限逐渐贴着负实轴从下方进入原点。第三步是求解穿越点。奈奎斯特图上最值得关注的穿越点有两类一是穿越虚轴的点令 (R(\omega)0)二是穿越实轴的点令 (I(\omega)0)。对上面的系统(R(\omega)) 恒为负不存在正实部区间(I(\omega)) 恒为负也没有实轴穿越点这说明曲线始终待在下半平面不会碰到 ((-1, j0)) 所在的区域稳定性判断一目了然。如果遇到更高阶的系统比如 (\frac{K}{s(s1)(s2)})求穿越点就需要解一个多项式方程。这种时候你可以借助数值求根工具也可以直接丢给后面的数值法。2.3 手工法的适用边界和常见错误手工解析法并不是用来求解复杂系统的它最适合传递函数阶次较低、结构标准的情况。它的最大价值在于建立直觉让你知道一条奈奎斯特曲线为什么长成那个形状起点终点方向怎么确定穿越点怎么求。我见过不少工程师直接用MATLAB画图画完却不敢确认结果对不对就是因为缺少这种直觉。反过来如果你能手工画出一阶惯性环节、积分环节、惯性加积分环节的草图再用软件验证工具出了问题你也能第一时间发现。手工法最常见的坑有三个一是不补负频率部分。奈奎斯特判据用的是从 (-\infty) 到 (\infty) 的完整曲线正频率部分画完还要知道关于实轴对称的负频率部分。二是忘了处理积分环节在原点附近的无穷大半圆。三是把 ((-1, j0)) 的包围方向搞反。建议你每次画完嘴上默念一遍“曲线从哪来、到哪去、兜了几个圈”比闷头算到底强得多。3. 方法二数值计算法——工程中最省心的高效路线3.1 用Python控制库五分钟出图只要是正经做控制的人电脑里大概率装了Python和control库。没有的话先装pip install control matplotlib然后几行代码就能出图以 (G(s)\frac{1}{s(s1)}) 为例import control as ct import matplotlib.pyplot as plt G ct.tf([1], [1, 1, 0]) # 分子1分母s^2s也就是s(s1) ct.nyquist_plot(G) plt.grid(True) plt.axhline(0, colorblack, linewidth0.8) plt.axvline(0, colorblack, linewidth0.8) plt.xlabel(Re) plt.ylabel(Im) plt.show()这里ct.tf([1], [1, 1, 0])创建传递函数分子系数是[1]分母系数是[1, 1, 0]对应 (s^2s)。ct.nyquist_plot会自动计算正的频率响应、负的频率响应并且在有积分环节时补好无穷远处的连接圆弧。老版本control库这个函数名是ct.nyquist如果你用的是老版本把函数名换成ct.nyquist即可。默认画出来的图包含正负两支频率响应曲线关于实轴对称看起来不像手工画的正频率单支但它对稳定性判断更严谨。3.2 想掌控每个细节就自己手动离散频率计算用现成库很省事但库帮你做了太多事反而容易让使用者失去判断。我建议你至少写一次手动版本import numpy as np import matplotlib.pyplot as plt w np.logspace(-3, 3, 2000) s 1j * w G 1 / (s**2 s) plt.plot(G.real, G.imag, labelpositive freq) plt.plot(G.real, -G.imag, labelnegative freq) plt.grid(True) plt.axhline(0, colorblack, linewidth0.8) plt.axvline(0, colorblack, linewidth0.8) plt.axis([-2, 1, -5, 5]) plt.show()这个手动版本的思路其实就是你在手工法里做的那些事把每个 (\omega) 代入 (G(j\omega))算出复数再一一点到复平面上。注意两点频率范围w我用的是从 (10^{-3}) 到 (10^3) 的2000个点原因是有积分环节时低频幅值很大如果从 (0) 开始会出现无穷大除零问题。绘图时用plt.axis限制了一下坐标范围不然起点附近的无穷大会把整个图压成一条线。我做了一个小反转正频率部分画完后负频率部分直接取共轭画在图上。理论上负频率就是正频率关于实轴的镜像这样画出来的结果和库函数基本一致。你可以用这个方法自查control库的默认输出。3.3 MATLAB一行命令与电子绘图中的隐藏参数MATLAB用户更幸福控制工具箱自带nyquist函数G tf(1, [1 1 0]); nyquist(G); grid on;这条命令画出的图默认包含正负频率和积分环节的无穷大半圆坐标轴范围也会自动匹配适合快速出正式报告图。如果你想拿数据做进一步分析可以用带返回值的调用方式G tf(1, [1 1 0]); [re, im, w] nyquist(G);这里返回的re和im是复数域上的频率响应实部虚部w是频率向量。需要注意MATLAB返回的re、im在SISO系统下是三维数组或带结构的类型直接画图前用squeeze或reshape处理一下。数值法看起来无脑实际坑点一点都不少。频率范围的选择直接影响曲线形状如果只选到 (\omega100)可能看不到高频终点完全回到原点如果低频只从 (\omega1) 开始含有积分环节的曲线起点就完全错过。我的建议是低频至少比系统最小转折频率小一个数量级高频至少比最大转折频率大一个数量级。系统模型变化大时就用np.logspace(-4, 4, 4000)这种更宽的扫描范围。3.4 数值法的几个高频翻车点第一个坑是采样点数量不够。奈奎斯特曲线在中频段往往弯得很厉害频率点太少曲线会变成折线穿越点也被藏掉。2000个点以上基本够用复杂系统建议上4000到10000个点。第二个坑是MIMO系统。直接用nyquist(G)画多输入多输出系统会输出多个子图每个输入输出通道各一套曲线图面非常拥挤。实际工程建议只关注关心的通道或者手动指定G[i, j]画某一对输入输出。第三个坑是不知道怎么验证数值结果。数值曲线看起来总是“很有道理”但它可能因为频率范围问题把关键的负频率部分掩盖了。我通常会在画完图后和手工法做个对照先确认几个特殊点起点在哪、终点趋向什么方向、是否穿越实轴。只有机器结果和手工直觉对上了我才放心用这张图做稳定性结论。4. 方法三由Bode图转换法——实测数据下最实用的替代方案4.1 为什么传递函数未知时第三条路反而香前面两种方法都建立在传递函数已知的基础上。但工程现场经常没这么理想你可能面对的是一个黑箱对象只有扫频仪测出来的Bode图数据也可能系统模型太复杂辨识出来一堆高阶项反而不如直接用实测频率响应。Bode图和奈奎斯特图本质是同一个复数频率响应的两种呈现方式。Bode图给出的是幅值 (|G(j\omega)|) 和相角 (\angle G(j\omega))而奈奎斯特图需要的实部和虚部可以用极坐标转直角坐标得到[ R(\omega) |G(j\omega)|\cos(\angle G(j\omega)) ][ I(\omega) |G(j\omega)|\sin(\angle G(j\omega)) ]如果你手头只有Bode图曲线没有现成数据文件依然可以人工取若干频率点读幅值读相位逐点转换成 (R) 和 (I)描点连线。这个方法看起来老派实际上在实测数据分析和现场校准时非常管用。4.2 Bode数据转奈奎斯特图的具体操作流程我按自己实际处理的流程拆解一下。第一步从Bode图上取点。取点的原则是低频段和相位快速变化段要密一些高频平坦段可以疏一点。工程上一般取10到20个点就够画趋势要精确定位穿越点就得在相角接近 (-180^\circ) 的区间多加几个点。第二步把幅值的dB值转成线性幅值。Bode图纵轴单位是dB奈奎斯特图的坐标是线性幅值转换公式[ A 10^{\frac{A_{dB}}{20}} ]第三步把相位从角度制转成弧度制计算实部虚部。可以在Excel里做也可以用Python写几行import numpy as np freq np.array([0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0]) # 单位可自定义 mag_db np.array([20.0, 17.0, 8.0, 0.0, -6.0, -14.0, -20.0]) phase_deg np.array([-95, -110, -140, -175, -195, -215, -230]) A 10 ** (mag_db / 20) phi np.deg2rad(phase_deg) Re A * np.cos(phi) Im A * np.sin(phi) import matplotlib.pyplot as plt plt.plot(Re, Im, o-) plt.grid(True) plt.show()这样得到的曲线和直接从传递函数计算的奈奎斯特图在工程精度范围内应该吻合。如果数据是从仪器里导出的只需要把读數值替换成CSV或Excel里的列。这里最容易出错的是相位符号。Bode图相位一般是负角度奈奎斯特图的虚部对应负相位时为负这个符号别搞反。4.3 不逐点转换也能用Bode特征快速估算稳定性逐点转换适合出图但现场判断稳定性时有一个更快的办法只关注相角等于 (-180^\circ) 的那个频率点通常记为相位穿越频率 (\omega_{pc})。在这个频率上Bode幅值是多少如果幅值大于 (0\text{ dB})意味着 (|G(j\omega)| 1)换算到奈奎斯特图上曲线在负实轴那一侧已经超出了以原点为中心、半径为1的单位圆也就是包住了 ((-1, j0)) 点附近。对于开环最小相位系统这种情况基本等价于闭环不稳定。举个例子某系统Bode图显示相位在 (\omega3\text{ rad/s}) 时穿越 (-180^\circ)此处幅值为 (6\text{ dB})那么奈奎斯特图上对应的点在负实轴上距离原点超过1((-1, j0)) 点被困在曲线左侧还是右侧你可以用笔算一下幅值 (10^{6/20} \approx 2)相位 (-180^\circ)该点坐标为 ((-2, 0))。从几何上看曲线在这个点附近走了一个大圈((-1, j0)) 点被包含进去的概率非常高。这正是Bode图能快速判断奈奎斯特图稳定性的底层逻辑。Bode转换法的局限也很明显实测数据的相角在高频段容易受噪声和采样率影响误差会被放大到实部虚部计算里非最小相位系统在相位穿越频率附近可能有不同走向不能直接套用最小相位的简单判断。所以这个方法适合工程初判正式结论还是要结合精确模型或时域响应验证。5. 三法对比、选型建议与避坑指南5.1 三种方法到底怎么选把三种方法的适用场景摆在一起看会更清楚方法适用场景精度对原理掌握要求工具依赖耗时手工解析法低阶系统、学习原理、快速草图中取决于特殊点选取高纸笔即可根据系统复杂度不定数值计算法传递函数已知的工程分析、正式报告高中Python/MATLAB短Bode图转换法传递函数未知、实测扫频数据、现场调试中高受取点影响中高Excel或脚本中我的实际建议是学习阶段三种都练至少把一个系统用三种方法各画一遍。工程上直接上数值计算法但出结论前用Bode特征点快速验证一下趋势。遇到黑箱系统老老实实走Bode转换法。三者不是互斥的而是互相校验的关系。5.2 频率范围和控制库的高级选项很多人在数值计算里只改传递函数不改频率范围这是不够的。Python control库画奈奎斯特图时可以通过传递omega参数指定频率范围w np.logspace(-2, 3, 5000) ct.nyquist_plot(G, omegaw)MATLAB对应写法w logspace(-2, 3, 5000); nyquist(G, w);选频率范围的原则我总结成一句话低频要低到能看到曲线起点高频要高到能看到曲线收进原点。对含有积分环节的系统低频再低也只是趋向无穷视觉上无所谓但稳定性判断时无穷大半圆的处理依赖工具补全所以你手动指定频率别从0开始避免除零。还有一个实用小技巧如果只想关心正频率半支可以在画图后只看实部虚部都较小的局部区域用坐标轴缩放观察穿越点附近。别想着让全局图把所有细节都展示出来奈奎斯特图通常要配合局部放大看关键段。5.3 稳定性判断口诀和自检流程根据我的经验初学者最容易死记公式然后套错符号。我整理了一个自检口诀用于每次画完图后检查自己的判断数一下P数一下圈顺时针绕一圈N加一算完Z看零没零不放心再用Bode裕度验一遍。具体操作流程可以这样确认开环右半平面极点数 (P)绝大多数工程对象 (P0)。看奈奎斯特曲线逆时针/顺时针绕 ((-1, j0)) 的净圈数确认正方向定义和你用的公式一致。代入 (ZPN)我这里用的定义是顺时针为正必须满足 (Z0)。用Bode图的幅值裕度和相位裕度交叉验证。相位裕度为正、幅值裕度大于1一般对应奈奎斯特图不包围临界点。这套流程我用了十几年能拦截掉绝大多数画图误判。特别是“工具画出来的图到底包没包住(-1,-j0)”这种问题用Bode裕度去对一下十有八九能发现频率范围设置不对导致的假包围。5.4 高频极点、非最小相位系统和离散系统最后提两个我踩过的进阶坑。第一个是开环含有右半平面零点的情况。非最小相位系统的奈奎斯特曲线走向往往和最小相位系统的直觉相反起点终点方向都不能直接套习惯。这种时候手工法的角色更重要你必须一步步推导实部虚部搞清楚相角变化趋势再让数值法给你出图。第二个是离散系统。离散传递函数画奈奎斯特图时频率变量不再是 (sj\omega)而是 (ze^{j\omega T})这时图形往往更复杂。好在control库和MATLAB都支持离散系统传入离散传递函数对象后nyquist_plot和nyquist都能正确处理。只是你看到图的坐标、频率范围含义和连续系统不完全一样别拿连续系统的经验硬套。说到底画奈奎斯特图的三种方法本质是同一个数学对象的三副眼镜。手工法让你看清楚原理数值法帮你高效出图Bode转换法让你在传递函数缺失时依然能干活。我个人的体会是你在实际项目中很少只用其中一种通常先用数值法出精确曲线再用手工法检查特殊点最后用Bode特征点给领导或者同事做个快速说明。这三种方法都熟练以后奈奎斯特图在你眼里就不再是一条需要背的曲线而是一套完整的稳定性分析思路。最后再分享一个小技巧每次画完图把鼠标或者光标移到曲线的起点、终点和所有实轴交点上看一眼这三个位置的坐标如果能跟你手算的数值对得上这张图基本不会错。希望这篇经验对你有所帮助。