MATLAB传递函数频域分析：Bode/Nyquist/Nichols与稳定裕度

简介面向控制系统与Simulink入门学习者的MATLAB小工具包聚焦于把S平面上的传递函数转换为直观的三维曲面图便于观察频率响应随频率与相位变化的特性辅助理解极点零点分布、稳定性与动态性能等抽象概念。压缩包共2个文件以m脚本和txt文本为主脚本承担传递函数定义、频率范围设置、响应计算与三维绘图等核心逻辑文本文件为使用许可说明整体约2KB轻量易取用。已有1822人学习说明该主题在课程实验与自学中具有一定关注度。读者可获得一份可直接参考的传递函数三维可视化实现用于对照tf函数建模、频响计算及surf、meshgrid绘图思路也可作为Simulink基础学习的补充素材适合需要把控制理论曲线与三维图形对应起来的中初级学习者。1. 三张图看清一个传递函数从 tf 对象到频域曲线某开发者把 Simulink 里的控制器参数改完打开 Bode 图一看幅频曲线往右挪了相位裕度却从 45 度掉到 12 度跟手算对不上。查了半天问题不在模型——他画图时手敲的分母系数漏了一项。在 MATLAB 里做传递函数绘图这类翻车太常见tf、zpk、ss 三种写法混着用Bode、Nyquist、Nichols 三张图各画各的中间没有一次交叉验证。这份资源围绕“一个传递函数、三张频域图”搭了一套 MATLAB 脚本与 Simulink 示例模型从建模、出图、算稳定裕度一路走到与 Simulink 线性化结果核对。刚接触 Simulink 基础、需要把频域指标落成具体数字的人能照着跑已经会画图、被单位、时延、高阶数值问题坑过的熟手也能在里面找到边界条件。2. 传递函数的三种写法与互转tf、zpk、ss 该用哪个传递函数在 MATLAB 里不是一段公式字符串而是一个带属性的对象。你敲tf([1 3],[1 4 5 0])存下来的是分子分母系数、输入输出名、时延、采样时间这一堆字段。后面画 Bode、算裕度、跟 Simulink 做线性化核对读的全是这些字段。所以第一步不是急着调bode而是先想清楚用哪种表示把这个对象建出来。选错了表示图能画出来但数字不一定可信。我见过把八阶系统的分母手写展开再喂给 tf 的图看着挺顺裕度算出来是负的仿真却稳得很——这种“图对、数错”的情况比直接报错难查得多。2.1 三种写法不是三个东西选型的三个判断点tf多项式、zpk零极点增益、ss状态空间描述的是同一个线性时不变系统差别在数值稳定性和信息密度不在功能。低阶、系数整齐、要给别人看分母长什么样用 tf 最直接。二阶以内的系统手写 tf 几乎不出错报告里贴出来也一目了然评审一眼能对上物理量。阶次过六或者零极点分布跨了好几个数量级我会强制换成 zpk。原因很实在tf 内部要靠多项式乘法展开分母展开本质是卷积系数一大一小挨着浮点有效位掉得很快。zpk 从头到尾只存零极点位置和增益不展开救的就是这一段精度。MIMO 系统、后面要接状态初值、要做状态反馈观测器用 ss。Simulink 的线性化输出天生就是 A、B、C、D 四个矩阵转成 tf 属于二次加工多一道工序就多一次失真机会。判断顺序我一般这么走先看输入输出是不是多于一对是就走 ss再看阶次是不是超过六是就走 zpk剩下的小系统才用 tf。这个顺序不是死规矩但能挡掉大部分精度问题。2.2 建立与互转一段能跑通的脚本包里第一节的脚本就是干这件事的建一个三阶系统派生三种表示再用频响数据验证互转有没有丢东西。别用肉眼判断“两条曲线看起来一样”要落到数字上。% 同一个三阶系统三种写法 num [1 3]; % 分子s 3 den [1 4 5 0]; % 分母s^3 4s^2 5s含一个积分环节 G_tf tf(num, den); % 多项式形式 G_zpk zpk(G_tf); % 零极点增益形式由 tf 直接派生 G_ss ss(G_tf); % 状态空间形式默认可控标准型 % 互转回来后逐频点比对不要靠眼睛 w logspace(-2, 2, 300); H_tf squeeze(freqresp(G_tf, w)); H_zpk squeeze(freqresp(tf(G_zpk), w)); H_ss squeeze(freqresp(tf(G_ss), w)); % 逐点相对误差用最大相对偏差衡量回转是否无损 err_zpk max(abs(H_tf - H_zpk) ./ max(abs(H_tf), eps)); err_ss max(abs(H_tf - H_ss) ./ max(abs(H_tf), eps)); fprintf(zpk 回转误差 %.3ess 回转误差 %.3e\n, err_zpk, err_ss);逻辑上先建 tf 作为基准再由它派生 zpk 和 ss保证三者描述的是同一个系统回转比对用freqresp在频域取复数响应逐点做相对误差。这里特意没用norm分母里有 s 因子系统在原点是极点H2 范数不存在直接norm(G_tf - tf(G_zpk))会抛错新手往往卡在这一步。参数上logspace(-2,2,300)覆盖 0.01 到 100 rad/s对一个三阶系统足够看清全部转折eps是防止某个频点响应恰好为零时除零freqresp返回的是 1×1×N 的三维数组squeeze压成一维才能跟频率向量对齐。误差量级正常在 1e-14 附近如果跑出来是 1e-3 甚至更大说明你的分母系数本身就有问题别急着往下画图。2.3 离散化与时延画图之前先变成有理函数连续系统的图好画麻烦的是离散控制器和纯滞后。这两类东西直接往 tf 里塞图能出来但读出来的裕度不可信。离散化先要定方法。c2d不写方法时默认零阶保持器如果 Simulink 模型里的离散模块用的是双线性两边模型根本不是同一个后面核对必然对不上。Ts 0.01; % 采样周期必须与 Simulink 里的 FixedStep 一致 Gd_zoh c2d(G_tf, Ts, zoh); % 零阶保持器最贴近实际数字控制器 Gd_tus c2d(G_tf, Ts, tustin); % 双线性变换低频段吻合好高频段会翘 Gd_mp c2d(G_tf, Ts, matched); % 零极点匹配只在通带内可信 % 时延不要直接把 exp(-sT) 塞进传递函数 G_delay G_tf; G_delay.InputDelay 0.05; % 50 ms 纯滞后 G_pade pade(G_delay, 4); % 4 阶 Pade 有理近似才能画 Bode 和 NyquistTs这个参数不是随便填的它决定了离散后频率轴能画到哪里上限是奈奎斯特频率pi/Ts图超过这条线就没有物理意义。pade的第二参数是近似阶数常用 3 到 5 阶阶数越高在低频段越准但会引入额外的右半平面零点高频段相位反而失真更厉害。判稳时只信 Pade 近似带宽以内的那一段超出的部分画出来好看不能拿来下结论。2.4 阶次上去以后poly 展开的精度账三个一阶环节串联手写分母看起来风平浪静实际系数动态范围能跨三个数量级。阶次再往上走问题就不是“图丑不丑”而是“分母还代不代表原系统”。s tf(s); G1 1/(s 1); G2 1/(s 10); G3 1/(s 100); G_conv G1 * G2 * G3; % 内部按多项式卷积展开分母 G_zpkm zpk(G1) * zpk(G2) * zpk(G3); % 零极点直接相乘不展开 [num_c, den_c] tfdata(G_conv, v); den_nz den_c(den_c ~ 0); spread max(abs(den_nz)) / min(abs(den_nz)); % 分母系数的动态范围 fprintf(分母系数动态范围 %.2e\n, spread);这段代码的重点是最后那个spread。它超过 1e6 就说明分母系数在双精度下已经开始互相“吃掉”有效位此时pole(G_conv)求出来的极点位置会有可见偏差再拿去画 Bode、算裕度误差会一路放大。tfdata(..., v)里的v是让返回值直接给数值向量而不是元胞数组省一次cell2mat。真遇到动态范围大的系统做法是全程用 zpk 建模、用 zpk 相乘只在导出报告给人看的时候临时转一次 tf。3. Bode 图落地频率范围、幅值单位与稳定裕度怎么设Bode 图是这三张图里最容易被当成“随便一画”的一张。默认频率范围常常把关键转折点挤在最左边一小段幅值单位默认是 dB 但有人拿线性值去读相位默认折回 ±180 度导致累计相位看不见。这三件事任意一件出错读出来的裕度就是错的。包里第三节的脚本把出图拆成四步先算频率范围再设绘图选项然后算裕度最后做多曲线对比和数据导出。顺序不能反频率范围定错了后面所有读数都建立在错误的坐标系上。3.1 频率范围不是拍脑袋由极点模值定上下限默认的bode(G)会自动选范围但对含积分环节的系统它经常从 1e-2 起步把低频段的 -20 dB/dec 斜率截掉一大块。手算过一次带宽之后我基本不再用默认范围。p pole(G_tf); p p(abs(p) 1e-9); % 剔除原点上的积分极点否则 min 会取到 0 if isempty(p) w_lo 1e-3; % 纯积分系统下限手工给 else w_lo 0.1 * min(abs(p)); % 低于最慢极点一个数量级 end w_hi 10 * max(abs(p)); % 高于最快极点一个数量级 w logspace(log10(w_lo), log10(w_hi), 500);逻辑是让频率窗口把全部极点模值包在中间两边各留一个数量级。低于最慢极点一个数量级能看到低频渐近线高于最快极点一个数量级能看到高频滚降的斜率对不对。abs(p) 1e-9这个过滤很关键积分环节的极点数值上是 0 或者极小的负数不剔掉会让log10(w_lo)直接变成负无穷。logspace的第三个参数是点数500 点对一张单系统 Bode 图足够密做数据导出时我会加到 2000 点避免 CSV 里曲线出现折线感。3.2 bodeoptions 里真正影响结论的四个参数bodeoptions有一堆属性绝大多数不影响你读出来的数字但这四个会。opts bodeoptions; opts.FreqUnits rad/s; % 与推导时用的单位一致别在 Hz 和 rad/s 之间来回换 opts.MagUnits dB; % 读裕度必须用 dB线性刻度容易看错量级 opts.PhaseUnits deg; opts.Grid on; opts.XLimMode manual; opts.XLim [w_lo w_hi]; opts.PhaseWrapping off; % 不折回 ±180看清累计相位 bodeplot(G_tf, opts);FreqUnits是最容易出事的一个。Simulink 里很多模块的频率参数按 Hz 填脚本里按 rad/s 推混一次就差 2π 倍穿越频率能差出一个数量级。PhaseWrapping设成off之后相位会一路往下走三阶系统能看到 -270 度这在判断“相位还差多少到 -180”时比折叠显示直观得多。XLimMode必须显式设成manual只给XLim而不改模式MATLAB 会按自动范围画你给的区间被忽略。这个坑我踩过一次调了半天以为选项名字写错了。3.3 用 margin 把曲线读成数字图是给人看的数字是给判据用的。margin一次给出四个量幅值裕度、相位裕度、相位穿越频率、幅值穿越频率。[Gm, Pm, Wcg, Wcp] margin(G_tf); Gm_dB 20*log10(Gm); % margin 返回的是线性幅值比报告里写 dB if isinf(Gm) fprintf(幅频不穿越 0 dB幅值裕度无穷大Wcp %.3f rad/s\n, Wcp); end T feedback(G_tf, 1); % 单位负反馈闭环 bw bandwidth(T); % 闭环 -3 dB 带宽 fprintf(Gm %.2f dB, Pm %.2f deg, Wcg %.3f, Wcp %.3f, BW %.3f\n, ... Gm_dB, Pm, Wcg, Wcp, bw);margin返回的Gm是线性比值不是 dB这一步换算漏掉的话报告里的幅值裕度会写成 3.16 而不是 10 dB差了一个对数关系。Wcg在一阶系统和过阻尼系统上常常是Inf或者NaN因为相位根本到不了 -180 度此时幅值裕度是无穷大属于正常现象不是脚本报错。bandwidth要喂闭环不是开环。开环的bandwidth在 MATLAB 里没有物理意义喂错了会得到一个跟闭环带宽差很远的数。3.4 多条曲线同图对比与数据导出调参阶段真正有用的不是单条曲线是同一坐标系下三四个增益版本的对比。G_list [G_tf, G_tf*1.2, G_tf*0.8]; % 标称、增益 20%、增益 -20% h bodeplot(G_list, opts); legend(标称, 增益 20%, 增益 -20%, Location, southwest); % 把幅频和相频导成表格方便贴进报告或做后处理 [mag, phase, wout] bode(G_tf, w); T_out table(wout(:), 20*log10(squeeze(mag)), squeeze(phase), ... VariableNames, {w_rad_s, mag_dB, phase_deg}); writetable(T_out, bode_export.csv);bode三输出的顺序是幅值、相位、频率跟bodeplot的句柄用法不一样别混。mag是线性幅值导数时手动20*log10一下导出的 CSV 才和图上看到的一致。squeeze在这里不能省1×1×500 的数组直接进 table 会变成一列元胞。参数常用取值改了会怎样FreqUnitsrad/s换成 Hz 后穿越频率数值差 2π 倍MagUnitsdB换成 abs 后 -3 dB 点要自己换算PhaseWrappingoff折回 ±180 会掩盖累计相位XLimModemanual不设 manualXLim 给的区间不生效4. Nyquist 与 Nichols包围圈判稳和等 M 圆读图Bode 图告诉你“离 -1 点还有多远”Nyquist 图告诉你“绕没绕过去”Nichols 图告诉你“闭环峰有多高”。三张图回答的是三个不同问题只画一张就下结论很容易在开环不稳定系统上翻车。包里第四节的脚本把三张图串起来跑同一组频率点、同一套单位画完互相核对。这一节的重点不是画图命令本身是判据怎么用、参数怎么设才能看清 -1 点附近那一小片区域。4.1 Nyquist 判稳Z P N符号和方向最容易记反Nyquist 判据的完整形式是Z P NZ 是闭环右半平面极点数P 是开环右半平面极点数N 是 Nyquist 曲线顺时针包围 -1 点的圈数。闭环稳定等价于 Z 等于零。麻烦在于 N 的符号约定。逆时针包围算正、顺时针算负不同教材写法不一样记串一次结论就反了。我的做法是不背符号直接用allmargin的Stable字段交叉确认。P sum(real(pole(G_tf)) 0); % 开环右半平面极点数 S allmargin(G_tf); % 结构体GainMargin、PhaseMargin、GMFrequency 等 fprintf(开环不稳定极点 %d 个各裕度对应闭环稳定标志\n, P); disp(S.Stable); % 手工取一段频率的实部虚部用于在图上定位 -1 点附近 [re, im] nyquist(G_tf, w); re squeeze(re); im squeeze(im); idx find(abs(re 1) 0.2 abs(im) 0.5); % 靠近 (-1,0) 的频点 fprintf(靠近 -1 点的频段%.3f ~ %.3f rad/s\n, min(w(idx)), max(w(idx)));allmargin的Stable是逻辑数组每个元素对应一组裕度下的闭环稳定性判断全为真才算稳。比手数圈数可靠得多尤其是曲线在 -1 点附近来回穿插的时候。最后那段找idx的代码是为了定位“该放大哪一段”下一小节会用到。4.2 nyquistoptions 与 -1 点附近的放大默认的 Nyquist 图会把正负频率两支都画出来镜像分支经常从 -1 点前面横穿过去把关键区域挡得严严实实。opts_n nyquistoptions; opts_n.FreqUnits rad/s; opts_n.ShowFullContour off; % 只画正频率段避免镜像分支遮挡 opts_n.XLimMode manual; opts_n.XLim [-2 1]; opts_n.YLimMode manual; opts_n.YLim [-1.5 1.5]; nyquistplot(G_tf, opts_n);ShowFullContour设成off只影响显示不影响判据计算这一点可以放心。XLim给[-2 1]、YLim给[-1.5 1.5]是把窗口压在 -1 点周围两倍范围曲线绕没绕过去一眼能看出来。范围开太大-1 点缩成一个像素等于白画。对含积分环节的系统Nyquist 曲线在低频端会趋向无穷远此时窗口要给得稍微宽一点否则曲线从图外飞进来看不出它从哪个方向绕。4.3 Nichols 图叠加 ngrid直接从图上读闭环谐振Nichols 图的横轴是开环相位、纵轴是开环幅值叠上ngrid之后曲线和等 M 圆的切点就是闭环谐振峰不用再算一遍闭环。nichols(G_tf); ngrid; % 等 M 圆 等相角线 hold on; nichols(G_tf * 1.2); % 增益抬 20% 后的曲线 hold off;读图方法曲线与 0 dB 横线的交点对应幅值穿越频率与 -180 度竖线的距离就是相位裕度曲线离哪条等 M 圆最近闭环谐振峰 Mr 就在那一档附近。增益抬高 20% 之后曲线整体上移如果它擦到了 3 dB 的等 M 圆说明闭环超调会明显变大。ngrid必须在nichols之后调用顺序反了网格会被曲线盖住。旧版本里ngrid会新开一个图窗新版直接叠加在当前坐标轴上跑出来不对先查版本。4.4 三张图交叉验证一张对照表同一组数据画三张图对不上的地方就是问题所在。包里给了一张对照表我自己也常用这套判据。关注量Bode 怎么看Nyquist 怎么看Nichols 怎么看幅值裕度相位 -180 处的负 dB 数曲线与负实轴交点到 -1 的距离曲线与 -180 竖线的纵向距离相位裕度0 dB 处的相位与 -180 的差单位圆交点角度曲线与 0 dB 横线的横向距离闭环谐振需另算闭环不直观与等 M 圆的最近距离稳定性两个裕度都为正是否顺时针包围 -1曲线是否越过 -180 且幅值大于 0 dB三张图读数如果互相矛盾九成是单位或者频率范围出了问题。先查FreqUnits再查MagUnits最后查是不是把开环和闭环搞混了。这个排查顺序比逐行看代码快得多。5. 避坑清单五条把图跑歪的血泪记录这一节列的五个问题全部是“脚本不报错、图也画出来了、结论是错的”。它们不会在命令行里给你红色警告只会安静地把一个错误的裕度写进报告。每条按现象、原因、解决三段记照着对一遍能省掉大半天。5.1 幅值单位写成 abs裕度差了一个对数关系现象bodeplot出来的曲线在 0 dB 附近穿过margin给的Gm却是 3.16报告里直接抄成“幅值裕度 3.16 dB”。原因MagUnits设成abs时图按线性刻度画而margin内部始终按线性幅值比计算两边一个用线性一个用对数读图的人把 3.16 当成 dB 了。解决绘图选项统一写opts.MagUnits dB代码里需要线性值时用10^(Gm_dB/20)换算别在两种量纲之间来回猜。所有打印语句里带单位的量变量名后缀都加上_dB。5.2 时延直接塞进 tf相位一路卷到 -720 度现象给系统加 50 ms 纯滞后Bode 相位从 -90 度一路掉到 -720 度Nyquist 图上出现好几圈螺旋判稳完全没法看。原因exp(-sT)不是有理函数画图时按无穷阶展开处理相位随频率线性下降时延越大绕的圈越多。Nyquist 判据里的“包围圈数”被时延产生的额外圈数污染数出来的 N 是假的。解决先pade(G_delay, 4)做有理近似再画图阶数取 3 到 5并且只在 Pade 近似有效的带宽内读裕度。判稳时如果时延占比很大老老实实回到时域仿真看阶跃响应别硬用频域判据。5.3 高阶系统用 tf 展开分母系数已经失真现象八阶系统Bode 图高频段幅值不降反升margin算出来的相位裕度是负数Simulink 仿真却稳得很。原因poly展开时系数动态范围跨了七八个数量级双精度有效位不够分母多项式本身已经不代表原系统pole求出来的极点位置都是错的。解决从建模开始就用 zpk环节相乘也用 zpk 相乘全程不展开只在导出报告给人看的时候临时转一次 tf。判断标准就是第 2 章那个spread超过 1e6 就别再用 tf 做计算。5.4 margin 的输入给成了闭环现象裕度算出来永远是无穷大或者相位裕度九十多度看着特别“稳”。原因margin要的是开环传递函数。把feedback(G,1)的结果喂进去环路已经闭合不存在“穿越 0 dB”这回事返回的Inf是数学上的必然不是系统真的好。解决margin(G_open)G_open是控制器乘被控对象闭环只用来算bandwidth、看阶跃响应、算灵敏度函数。脚本里给变量起名时就把_open和_cl分开别用G1、G2这种名字。5.5 离散系统没指定 c2d 方法图和仿真对不上现象脚本画的 Bode 和 Simulink 里 Linear Analysis 出的曲线在 0.5 倍奈奎斯特频率附近差了 3 dB 以上低频段却几乎重合。原因c2d不写方法时默认zoh而 Simulink 里的离散模块用的是双线性或者实际控制器就是双线性离散两边根本不是同一个模型。低频段两种方法吻合高频段才分家所以容易漏看。解决c2d显式写方法和 Simulink 模型里的离散模块设置逐项对齐核对时两边用同一个Ts。图只看奈奎斯特频率以内的部分超过pi/Ts的曲线是数学延拓没有物理意义。6. 从脚本到 Simulink 联调用 linmod 校核三条曲线脚本里的传递函数是手写的Simulink 模型是搭出来的两者是不是同一个东西得用数据说话。我一般用linmod在工作点附近把模型线性化成状态空间再和脚本里的 tf 在同一个频率向量上比频响。load_system(plant_demo); % 打开示例模型不弹窗口 sys_lin linmod(plant_demo); % 工作点附近线性化返回 a/b/c/d 结构体 G_lin ss(sys_lin.a, sys_lin.b, sys_lin.c, sys_lin.d); w logspace(-2, 2, 300); H_script squeeze(freqresp(G_tf, w)); H_lin squeeze(freqresp(G_lin, w)); dev_dB 20*log10(abs(H_lin - H_script) ./ max(abs(H_script), eps)); [dev_max, i_max] max(dev_dB); fprintf(最大偏差 %.2f dB %.2f rad/s\n, dev_max, w(i_max));linmod只对连续模型有效模型里只要有离散模块就会报错或者给出错误结果这种情况要换成linearize配合linio指定分析点。偏差在 -40 dB 以下说明两边模型一致偏差集中在某个频段多半是某个环节的饱和或者限幅在工作点附近被激活了线性化结果不再代表原系统。调参阶段我还会把裕度扫成一张表增益从 0.5 倍扫到 2 倍每步算一次margin把Pm和bw并排列出来比一张一张图翻快得多。从那以后我每次出频域图之前都强制走一遍三张图核对Bode 看形状和裕度Nyquist 看包围方向Nichols 看闭环峰三个对不上就先查单位和开环闭环。这套习惯救过我至少两次一次是分母漏项一次是c2d方法没对齐。希望帮到你。本文还有配套的精品资源点击获取