做控制系统分析的人十有八九都经历过这种时刻拿到一个开环传递函数第一反应是打开Bode图想看幅频特性和相频特性得来回切换两张图要估个稳定裕度还得对着坐标轴反复确认。后来习惯了奈奎斯特Nyquist图才发现这家伙一张图就能把频率从0到∞的幅值与相位变化全部装进去而且G(jω)曲线和(-1, j0)点的位置关系直接对应闭环稳定性在很多工程场景下比Bode图更直觉。这篇东西想做的事很具体用同一个开环系统从零起步演示三种画奈奎斯特曲线的方法。第一种是最老派但极其有用的一招——手工推导把实部虚部分离出来沿着ω从0到∞的方向把曲线骨架画出来第二种是MATLAB的nyquist命令适合快速出图确认结果第三种是用python-control库适合想把这步绘图塞进自动化脚本或科研流水线的人。三种方法的步骤、原理和坑我都会写出来不管你是刚学自动控制原理的学生还是常年跟传递函数打交道的工程师应该都能找到合适的一种。1. 一张奈奎斯特图把开环频率特性从0到∞全讲清楚1.1 为什么有了Bode图还要用它奈奎斯特图的不可替代之处Bode图的核心优势是把幅值和相位拆成两张图横轴是频率刻度方便做渐近线近似所以大多数教材都先讲Bode。但它有个绕不开的短板幅值特性和相位特性分开后你必须同时盯着两张图才能判断系统稳不稳而且相位穿越-180°的位置、幅值穿越0dB的位置都要靠眼睛在两图之间“对坐标”。奈奎斯特图换个思路把频率ω当作隐变量横轴画Re[G(jω)]纵轴画Im[G(jω)]一条曲线同时包含幅值和相位信息。ω从0增加到∞曲线逐渐扫过复平面幅值就是曲线上点到原点的距离相位就是该点与原点连线的夹角。判断闭环稳定性时只需要看曲线绕没绕(-1, j0)点这就是幅角原理在工程上的落地。很多人第一次看奈奎斯特图会觉得抽象因为它把一个频率轴折叠成了曲线上的“进度条”。但换个角度想Bode图是“幅值和相位分别对着频率展开”奈奎斯特图是“把频率当作参数让复数轨迹自己走一遍”信息量其实一模一样只是投影方式不同。等你习惯了用“起点、终点、穿越点”三个要素去读图大部分系统都能在几秒内看出稳定裕度的量级。1.2 贯穿全文的演示系统为什么选这个传递函数全文统一使用这个开环传递函数G(s) 2 / [s(s1)(0.5s1)]这不是随便选的。它有三个典型特征非常适合用来演示三种绘制方法包含一个积分环节分母里的sω→0时幅值趋向无穷大是标准I型系统能逼你处理“曲线从无穷远出发”的情况。两个惯性环节的时间常数分别是1和0.5频率范围跨越两个数量级画图时能明显看到曲线转折。曲线会在某个频率穿越实轴穿越点和(-1, j0)的相对位置直接对应稳定裕度方便三种方法之间互相验证。先说结论后面三种方法都会回到这个结果穿越实轴的频率ωc≈1.414 rad/s对应实部坐标Re(ωc)-2/3。因为增益是2距离临界点(-1, j0)还差一段系统稳定如果把增益放大到3曲线会恰好穿过-1点系统进入临界稳定状态。2. 方法一手工推导慢慢画把曲线骨架算明白2.1 把s换成jω再分离实部与虚部手工绘制看起来繁琐但它能让你真正理解曲线形态的来源尤其是当工具给出一个奇怪图形时只有手算过的人才能判断是系统本身特性还是参数错误。第一步令sjω代入开环传递函数G(jω) 2 / [jω(1jω)(10.5jω)]接下来需要展开分母。我见过很多同学在这一步直接用数值代其实先做代数化简更省力。先把两个一阶因子相乘(1jω)(10.5jω) 1 j1.5ω - 0.5ω²再乘以外面的jωjω(1 j1.5ω - 0.5ω²) -1.5ω² j(ω - 0.5ω³)所以G(jω) 2 / [-1.5ω² j(ω - 0.5ω³)]要分离实部虚部就得消掉分母里的虚数。把分子分母同时乘以分母的共轭这是复数运算的标准操作类似把分母有理化G(jω) 2[-1.5ω² - j(ω - 0.5ω³)] / [(-1.5ω²)² (ω - 0.5ω³)²]分母展开后为(-1.5ω²)² 2.25ω⁴(ω - 0.5ω³)² ω² - ω⁴ 0.25ω⁶两项相加得到ω² 1.25ω⁴ 0.25ω⁶。于是实部虚部分别为Re -3ω² / (ω² 1.25ω⁴ 0.25ω⁶) -3 / (1 1.25ω² 0.25ω⁴)Im -2(ω - 0.5ω³) / (ω² 1.25ω⁴ 0.25ω⁶) -2(1 - 0.5ω²) / [ω(1 1.25ω² 0.25ω⁴)]到这一步后面的工作就变成纯代数和函数分析了。注意Re的表达式是“负”的说明曲线主体位于左半平面Im在ω√2时为负ω√2时为正说明曲线会从第三象限穿到第二象限。2.2 三个关键点起点、终点和穿越实轴的位置手绘不用描几百个点抓住几个关键点就能定出曲线骨架。第一个关键点是ω→0的起点。看Re和Im的表达式当ω→0时Re→-3Im的分子趋近-2而分母中ω的一次项趋于0所以Im→-∞。这意味着曲线从(-3, -∞)出发也就是第三象限的“最深处”。这个“实部有限、虚部负无穷”的组合正是I型系统的标志。第二个关键点是ω→∞的终点。当ω很大时Re≈-3/(0.25ω⁴)→0⁻Im≈-2(-0.5ω²)/[ω·0.25ω⁴]4/ω³→0⁺。所以曲线最终从第二象限方向收敛到原点也就是(0⁻, 0⁺)。第三个关键点是曲线与实轴的交点。令Im0得到1-0.5ωc²0即ωc√2≈1.414 rad/s。代回Re的表达式Re(ωc) -3 / (1 1.25×2 0.25×4) -3 / 4.5 -2/3所以穿越点在(-2/3, 0)。这个点非常重要它决定了系统离(-1, j0)还有多远直接对应增益裕度。2.3 选几个中间频率描点把曲线趋势画出来有了三个关键点再补几个中间频率的点曲线形态就出来了。我算了几组数据ω (rad/s)ReIm0-3-∞0.5-2.26-2.641-1.2-0.41.414-0.66702-0.30.14-0.0350.041∞0⁻0⁺把这几个点按ω增大的顺序连起来从(-3, -∞)出发先向右上方走经过(-2.26, -2.64)再经过(-1.2, -0.4)此时曲线还在第三象限但已经明显抬头到达(-2/3, 0)穿越实轴进入第二象限然后继续向右逼近原点在原点处收尾。整体形状像一条“勺子”的侧面。实际手绘时低频段因为Im值变化剧烈可以多取两个点高频段Re和Im都快速趋近0只需确认趋势即可。画完以后在曲线上标注ω增加方向的箭头把ω0、ω∞、穿越点三个位置标清楚这张图就具备工程可读性了。2.4 手绘结果直接读出稳定裕度手工画完曲线立刻能干一件事判断闭环稳定性。开环系统极点分别是0、-1、-2都在虚轴上或左半平面所以开环P0。按Nyquist判据如果从ω0到∞的曲线不包围(-1, j0)点闭环稳定如果绕过去就存在右半平面闭环极点。我们画的曲线穿越实轴于-2/3在-1点的右侧没有包围临界点所以K2时闭环稳定。如果增益从2增大到K穿越点的实部坐标变成-K/3推导时把分子2换成K即可临界条件为-K/3-1即K3。于是立刻得到一个重要参数增益裕度Gm3/21.5约为3.52dB。这就是手绘的价值。很多教材只教你手绘不解释这步能读出什么等你真正需要即时判断一个传递函数稳不稳的时候会发现这张草图画起来只需要五六分钟却比翻箱倒柜找仿真软件快得多。3. 方法二MATLAB一条命令出图但别把判断交给默认设置3.1 用tf对象建模然后nyquist一行绘图MATLAB是控制工程里最常见的工作环境画Nyquist图的命令也非常短s tf(s); G 2/(s*(s1)*(0.5*s1)); nyquist(G); grid on;执行之后MATLAB会弹出一张Nyquist图横轴实部、纵轴虚部自动选择频率范围。对于简单系统这个图基本能看但我会直接说默认设置并不总是够用尤其是面对带积分环节的系统时低频段经常被截断曲线看起来就像“凭空出现在半路”很容易误判。还有一个细节MATLAB的nyquist命令默认会同时画出ω从-∞到0的对称曲线和ω从0到∞的主曲线频率增加的走向用箭头标出。图形界面上可以用Data Cursor在曲线上取点能直接读到某个频率下的实部和虚部验证手工计算结果很方便。3.2 如何挑选频率范围和网格密度按下图我建议立即放弃默认频率范围改用对数间隔手动指定w logspace(-2, 2, 800); nyquist(G, w); xlim([-5 1]); ylim([-4 1]); axis equal; grid on;为什么用logspace奈奎斯特曲线在低频和高频段的形态变化差异很大线性间隔会导到高频点数太稀疏。logspace(-2,2,800)意味着从0.01到100 rad/s均匀分布在对数坐标上800个点对于这种三阶系统足够密曲线插值平滑也能避免穿越点附近因采样稀疏而错过。设置axis equal这一点经常被忽略。如果横纵轴比例不一致曲线的视觉形态会失真穿越点的相对位置也会产生错觉。工程图虽然不需要像几何制图那样严格但为了判断-1点到底在曲线内侧还是外侧保持等比例是最安全的做法。3.3 从图和数据里精确读取穿越频率看曲线大致够用但如果要写进报告最好还是从数值上提取。nyquist命令其实会返回频率响应数据[re, im, wout] nyquist(G, w); re squeeze(re); im squeeze(im); d diff(sign(im)); idx find(d ~ 0); for k 1:length(idx) i idx(k); wc(k) 0.5*(wout(i) wout(i1)); ReCross(k) 0.5*(re(i) re(i1)); end这段逻辑很简单虚部符号变化的位置就是曲线穿越实轴的位置。之所以用符号变化而不用abs(im)最小值是因为符号变化更鲁棒不会把“接近实轴但没有穿越”的点误判成穿越点。对演示系统跑这段代码得到wc≈1.414 rad/sReCross≈-0.667和手工推导完全一致。当然如果只是想快速验证直接在图上用Data Cursor点穿越位置也可以但自动化提取的价值在于当你需要对几十个传递函数批量做判断时手点是不现实的。3.4 结合margin函数交叉验证稳定裕度MATLAB还提供margin函数直接计算稳定裕度[Gm, Pm, wcg, wcp] margin(G);margin返回线性增益裕度Gm、相位裕度Pm、穿越频率wcg和截止频率wcp。对演示系统Gm≈1.5wcg≈1.414 rad/s正好和Nyquist曲线穿越点对应的临界增益3相符。这里建议养成交叉验证的习惯画完nyquist后顺手跑一下margin两者对不上说明哪里有误解对上了就相当于多了一层保险。有一点要注意margin对积分环节和阻尼很小的系统有时会警告或者给出奇怪结果这时候回到Nyquist图手工看穿越点反而是最可靠的兜底方案。4. 方法三python-control开源绘制把奈奎斯特图纳入自动化脚本4.1 安装与建模不同版本的API差异先确认如果你不想依赖商业软件或者希望把绘图和分析写进Python脚本python-control库是目前最顺手的开源方案。安装很简单pip install control绘图本身依赖numpy和matplotlib如果还想用更多的频域分析函数可以额外安装slycot但画Nyquist图和计算频率响应不需要它。python-control的API在不同版本之间有变化这是我特别想提醒的。早期版本用ct.nyquist(G)新版本有些改成了ct.nyquist_response和ct.nyquist_plot。写脚本之前先看一眼版本和帮助文档import control as ct print(ct.__version__) help(ct.nyquist_response)实际建模代码和MATLAB很像import control as ct s ct.TransferFunction.s G 2 / (s * (s 1) * (0.5*s 1))这里用G 对象就建好了连续时间传递函数。注意分母展开要自己写python-control的TranferFunction构造函数也接受多项式系数但用(s1)这种因式相乘的写法更直观也更容易检查错误。4.2 核心绘制命令与数据返回新版推荐的绘图方式是import numpy as np import matplotlib.pyplot as plt omega np.logspace(-2, 2, 1000) re, im, _ ct.nyquist_response(G, omega)nyquist_response返回实部、虚部数组和频率数组直接拿去做自定义绘图自由度比MATLAB高很多plt.figure(figsize(6, 6)) plt.plot(re, im, b-, linewidth1.5, labelω: 0→∞) plt.plot(re, -im, b--, linewidth0.8, labelω: -∞→0负频部分) plt.plot(0, 0, ko, label原点) plt.plot(-1, 0, r, markersize12, label(-1, j0)) plt.axhline(0, colorgray, linewidth0.5) plt.axvline(0, colorgray, linewidth0.5) plt.xlim([-5, 1]) plt.ylim([-4, 1]) plt.gca().set_aspect(equal) plt.grid(True) plt.legend() plt.show()负频率部分的虚线是根据对称性手动补的对实系数传递函数默认会画出正负频率两侧但有时你看不到负频曲线所以在这个演示里我用plt.plot(re, -im)直接画对称曲线保证视觉效果完整。如果你只有老版本也可以完全不用nyquist_response直接从传递函数计算频率响应omega np.logspace(-2, 2, 1000) H np.squeeze(G(1j*omega)) re H.real im H.imag这行代码本质是“让计算机代做我们手算的G(jω)”逻辑最透明出错也容易排查。4.3 从频率响应数据里自己找穿越点不依赖额外工具包时穿越点计算可以用线性插值法实现sign_change np.where(np.diff(np.sign(im)) ! 0)[0] for i in sign_change: w1, w2 omega[i], omega[i1] im1, im2 im[i], im[i1] re1, re2 re[i], re[i1] t im1 / (im1 - im2) wc w1 t*(w2 - w1) rec re1 t*(re2 - re1) print(f穿越点: Re{rec:.4f}, ω{wc:.4f} rad/s)核心思路是假设频率点之间实部虚部近似线性变化然后按虚部过零的比例插值。对演示系统跑出来Re≈-0.6666ω≈1.414基本就是解析值。这个方法的局限在于频率点间隔必须足够密否则插值误差会变大但1000个对数点对这种低阶系统绰绰有余。4.4 与Bode图叠加形成一条流水线python-control的另一个优势是能方便地把Bode图和Nyquist图拼在一起并且把稳定裕度计算整合进同一套脚本。比如我在日常分析中经常写这样的流程G 2 / (s * (s 1) * (0.5*s 1)) fig, (ax1, ax2) plt.subplots(1, 2, figsize(12, 5)) ct.bode_plot(G, omega, axax1) ax1.set_title(Bode Plot) re, im, _ ct.nyquist_response(G, omega) ax2.plot(re, im, b-, linewidth1.5) ax2.plot(re, -im, b--, linewidth0.8) ax2.plot(-1, 0, r, markersize12) ax2.axhline(0, colorgray, linewidth0.5) ax2.axvline(0, colorgray, linewidth0.5) ax2.set_aspect(equal) ax2.set_xlim([-5, 1]) ax2.set_ylim([-4, 1]) ax2.grid(True) ax2.set_title(Nyquist Plot) plt.tight_layout() plt.show()这套流程放进Jupyter Notebook之后所有分析都有原始代码可以追溯改个参数重跑一遍就行特别适合参数扫描和批量对比。5. 三种方法放在一起怎么选5.1 横向对比速度、精度、适用场景三种方法没有绝对的“谁更好”它们适配的场景不同。我列了一个对比表对比项手工推导MATLABpython-control上手门槛高需要复数运算功底低一条命令中等需要Python基础时间成本5-15分钟几秒几秒精度取决于关键点选取高高批量自动化不适合可脚本化很方便对原理理解帮助最大一般一般授权成本无商业授权开源免费典型场景教学、考场、快速判断工程报告、学术论文自动化脚本、科研复现如果你还在学自动控制原理我的建议很直接先把方法一练熟至少练五个不同类型的系统再开始用工具。否则工具画出来的曲线对你来说就是一个没有依据的图形遇到异常也无从判断。如果你是在企业做工程确认或出报告MATLAB是主流选择因为格式标准、资料多、团队协作方便。如果你是做科研或者喜欢把一切都脚本化python-control是不二之选。它最大的优势是能和机器学习、数据分析的其他Python包写进同一个流程自由度极高。5.2 手工绘图仍是数字工具的地基这不是客套话。我见过很多能熟练使用MATLAB的工程师但在评审会上被问“为什么曲线是这个形状”时回答不出来。根本原因就是从没手算过不知道曲线形态是由哪些极点和零点决定的。手工绘图的价值在于建立“因果链”看到分母有积分环节就知道曲线低频段趋向无穷看到一对惯性环节可以画到什么位置就能预估穿越点大概在哪个象限一旦工具给出的曲线形态和预判不符立刻意识到模型可能设定有误。这种“预判-验证”的能力只有亲手推导过实部虚部才养得出来。5.3 推荐一套从手算到工具的验证工作流下面是我个人比较常用的流程分享出来给你们参考拿到传递函数后先不碰计算机手算Re、Im表达式。算三个关键量ω→0的起点、ω→∞的终点、Im0的穿越频率和穿越点实部。用这些量大致画出曲线骨架标注(-1, j0)点位置。再用MATLAB或Python画精确曲线检查默认图是否和手绘骨架一致。提取穿越频率和手算值对比。用margin或自编的插值代码算增益裕度。如果存在矛盾回头检查是代数计算出错还是工具频率范围截断导致图形不完整。这套流程看起来多了一步“手算”但它能帮你提前发现问题避免拿到一张看似精美但方向错误的图。6. 实操中绕不开的细节和坑6.1 频率范围没选对图全部挤成一团这是最常见的问题。自动选择的频率范围往往以中频为主低频段不够低、高频段不够高。对于演示系统如果默认频率从0.1开始低频部分的“尾巴”就会被截掉曲线起点看起来在某个有限位置很容易误以为系统是0型系统。我建议绘制时把频率范围至少覆盖到穿越频率前后各两个数量级。演示系统的穿越频率约1.4 rad/s用logspace(-2, 2)就留出了0.01到100 rad/s的余量。如果系统有更大的时间常数低频端点还需要再往左扩。6.2 箭头和ω方向画出来不代表读对奈奎斯特曲线是有方向的曲线方向代表频率增大的方向。手绘时一定要在曲线上标箭头工具绘图默认也会标。读图判断稳定裕度时方向决定了你判断“正穿越”还是“负穿越”。有次我在评审PPT上看到一张没有箭头的Nyquist图提问者问“曲线从哪边绕过去的”全场沉默。这种细节在自动控制这门课里不是小事方向缺失会让判据变得不可用。自己用matplotlib画图时如果nyquist_response返回的数据里没带箭头可以用plt.quiver在几个点手动加箭头标注ω0.5、1、1.414、2这些关键频率位置。6.3 有积分环节时别忘了补那笔“半圆弧”很多初学者画完ω从0到∞的曲线就完事了直接在图上判断稳定性。但对于含有积分环节的I型、II型系统按Nyquist判据的标准路径ω从-∞到0、0到∞之间需要补一个半径无穷大的一段圆弧这段圆弧的走向和角度取决于系统型别。这段补线在图上通常画成从实轴某点出发的大圆弧。它的作用是把曲线在原点附近断开的部分连接起来让包围次数的计算有闭环逻辑。实际判断稳定性时由于补线位于无穷大半径处大多数情况下绕不到-1点但如果系统型别很高比如II型以上补线方向稍有偏差就会影响包围次数。我在手工画图时会在草图右方用虚线标注“此处补无穷大半径圆弧”提醒自己判据的完整性用工具画图时则单独把ω从0到∞的主曲线提取出来再按标准规则把补线手动画上去。6.4 从图中读数时容易忽视的坐标比例问题如果横纵轴比例不同曲线的视觉形态会变形最直接的影响是低估或高估穿越点到-1点的距离。MATLAB默认不保证等比例python-control的matplotlib默认更不会。所以我在前面反复提到axis equal或set_aspect(equal)。设置等比例之后你会看到一个直观效果演示系统的穿越点-2/3到-1点的距离是一段约1/3的线段这个视觉比例和实际增益裕度1.5是对应的。如果纵轴被压扁这1/3的间距可能看起来只有一丁点容易让人产生“快不稳定了”的错觉反过来如果横轴被拉伸又会觉得裕度很大。6.5 关于插值和采样的一个提醒提取穿越点时不要用单个最小值点也不要用太粗的频率采样。前者会把负频段或接近实轴但没有穿越的曲线误判为穿越点后者会让插值出来的穿越频率误差偏大。对数间隔、1000个点对绝大多数低阶系统都是稳妥的组合。如果遇到高阶系统或曲线拐弯很急的系统可以先把数据画出来观察穿越点附近有没有明显曲率变化再把该局部的频率点加密重跑一遍。这套“粗扫-加密-提取”的思路和实验测试里的两遍测量法是一个道理。我自己的操作习惯是接到一个不熟悉的传递函数先在纸上花几分钟手算关键点再用python-control绘制完整曲线做复核。当手算的穿越点和插值提取的穿越点对上时这个系统的脾气基本就摸清了。这个方法不算快但能让人对每条画出来的曲线都有底气而不是对着屏幕上一堆像素半信半疑。