别再死记硬背了！用MATLAB/Simulink手把手教你画Bode图和Nyquist曲线（附代码）

工程实践指南：用MATLAB/Simulink玩转频率特性分析与系统稳定性验证

在自动控制系统的设计与调试过程中，频率特性分析是不可或缺的核心技能。无论是评估系统稳定性还是优化动态性能，Bode图和Nyquist曲线都为我们提供了直观的工程语言。但对于初学者而言，从理论公式到工程实践往往存在一道难以跨越的鸿沟——如何将课本上的传递函数转化为可操作的工程工具？这正是MATLAB/Simulink大显身手的舞台。

传统教学中，学生需要花费大量时间手工绘制渐近线、计算相位裕度，这种重复性劳动不仅效率低下，更可能掩盖了对核心概念的理解。现代工程软件的价值在于：它能将我们从繁琐的计算中解放出来，专注于系统行为的本质分析。本文将带您跨越理论与实践的断层，通过具体案例演示如何用MATLAB脚本和Simulink模型实现从传递函数到稳定性分析的全流程自动化。

1. 频率特性分析的工程意义与MATLAB实现基础

频率特性分析之所以成为控制工程师的必备技能，源于其独特的工程价值。与时间域分析相比，频率响应法不需要求解微分方程，仅通过正弦信号激励就能全面评估系统性能。这种"以简驭繁"的特性使其在工业现场大放异彩——我们无需知道系统的精确数学模型，通过实验测量频率响应就能进行控制器设计。

MATLAB环境为这种分析提供了完美的工作平台。其控制系统工具箱(Control System Toolbox)包含专为频率响应分析优化的函数集，从基本的bode()、nyquist()到高级的margin()、allmargin()，形成了完整的分析链条。考虑一个典型二阶系统：

% 定义二阶系统传递函数 wn = 10; % 自然频率(rad/s) zeta = 0.5; % 阻尼比 G = tf(wn^2, [1 2*zeta*wn wn^2]); % 生成Bode图 figure bode(G) grid on title('二阶系统Bode图')

这段简洁的代码背后，MATLAB自动完成了从频率点采样到幅值/相位计算的复杂过程。工程师可以直观观察到谐振峰值、带宽等关键特征，而无需关心底层计算细节。对于实际工程系统，这种自动化分析的优势更为明显——当系统阶数超过四阶时，手工计算几乎不可行，而MATLAB依然能提供精确的频率特性曲线。

频率特性曲线的工程解读要点：

• 增益裕度：系统在变为不稳定前可增加的增益量，反映绝对稳定性
• 相位裕度：系统在失稳前可容忍的相位滞后，反映相对稳定性
• 带宽频率：幅值衰减至-3dB时的频率，表征系统响应速度
• 谐振峰值：幅频特性最大值，指示系统振荡倾向

2. 从传递函数到Bode图的完整生成流程

Bode图作为频率特性最常用的可视化工具，其绘制过程在MATLAB中已被极大简化。但要想获得具有工程价值的分析结果，仍需遵循系统化的操作流程。我们以一个工业温度控制系统为例，演示专业级的Bode图生成方法。

假设系统开环传递函数为： $$ G(s) = \frac{5(0.5s+1)}{s(2s+1)(0.1s+1)(0.02s+1)} $$

对应的MATLAB实现为：

% 构建传递函数模型 num = 5*[0.5 1]; den = conv([1 0], conv([2 1], conv([0.1 1], [0.02 1]))); G = tf(num, den); % 定制化Bode图绘制 figure opts = bodeoptions; opts.FreqUnits = 'Hz'; % 频率单位设为Hz opts.MagUnits = 'dB'; % 幅值单位设为dB opts.PhaseWrapping = 'on'; % 相位自动解缠绕 opts.Grid = 'on'; [mag,phase,w] = bode(G, opts); % 计算稳定裕度 [Gm,Pm,Wcg,Wcp] = margin(G); disp(['增益裕度: ', num2str(20*log10(Gm)), ' dB at ', num2str(Wcg), ' rad/s']) disp(['相位裕度: ', num2str(Pm), ' deg at ', num2str(Wcp), ' rad/s']) % 标注关键参数 subplot(2,1,1) title(['Bode图 - 相位裕度: ', num2str(Pm), '°']) hold on semilogx([Wcp Wcp], [-100 0], 'r--') text(Wcp, -3, ['\omega_c = ', num2str(Wcp), ' rad/s'])

专业级Bode图绘制技巧：

对于包含多个子系统的复杂模型，MATLAB支持更灵活的分析方式。例如比较控制器参数变化对系统的影响：

% 定义不同PID参数下的系统 Kp_values = [1 2 5]; figure hold on for Kp = Kp_values C = pid(Kp, 0.1, 0.01); bode(C*G) end legend(strcat('Kp=', string(num2cell(Kp_values)))) title('不同比例增益下的Bode图对比')

这种参数敏感性分析在设计阶段极为宝贵，它能帮助工程师快速评估不同控制策略的效果，而传统手工绘图方法几乎无法实现这种多维比较。

3. Nyquist曲线绘制与稳定性判据实践

Nyquist稳定性判据是频率域分析中最强有力的工具之一，它通过开环频率特性判断闭环系统稳定性。MATLAB的nyquist()函数虽然能自动生成曲线，但要准确解读结果仍需掌握特定技巧。

考虑一个条件稳定系统： $$ G(s) = \frac{10(s+1)}{s(s-1)(0.1s+1)} $$

其Nyquist分析代码如下：

% 条件稳定系统示例 G = tf(10*[1 1], conv([1 -1 0], [0.1 1])); % 基本Nyquist图 figure nyquist(G) axis([-15 5 -10 10]) % 手动调整坐标范围 grid on % 增强型Nyquist分析 w = logspace(-2,3,1000); [re,im] = nyquist(G, w); re = squeeze(re); im = squeeze(im); figure plot(re, im) hold on plot(-1, 0, 'ro', 'MarkerSize', 8) % 标记临界点(-1,j0) xlabel('实部') ylabel('虚部') title('增强型Nyquist图') grid on % 自动稳定性判定 [Gm,Pm,Wcg,Wcp] = margin(G); if Pm > 0 && Gm > 1 disp('系统稳定') else disp('系统不稳定') end

Nyquist判据工程实践要点：

1. 曲线走向：频率增加方向必须明确标注，通常用箭头表示
2. 临界点识别：准确标记(-1,j0)点位置，这是稳定性判断基准
3. 奇异点处理：对于开环不稳定系统，需补足无穷大半圆路径
4. 穿越计数：计算曲线围绕(-1,j0)点的净环绕次数

当处理非最小相位系统时，Nyquist分析尤为有用。例如某飞行控制系统开环传递函数包含右半平面极点：

% 非最小相位系统示例 G = tf([1 -2], conv([1 1 0], [0.5 1])); % 高级Nyquist分析 figure nyquist(G) hold on % 绘制单位圆参考线 theta = linspace(0, 2*pi, 100); plot(-1 + cos(theta), sin(theta), 'k--') title('非最小相位系统Nyquist分析')

这种情况下，传统Bode图难以直接判断稳定性，而Nyquist曲线能清晰显示环绕情况。工程师可以直观看到曲线是否包围临界点，以及包围的方向和次数，这是理论教学与工程实践完美结合的典范。

4. Simulink中的频率特性实时验证方法

虽然MATLAB脚本提供了精确的分析工具，但在实际工程中，我们常需要在Simulink环境下进行系统级验证。Simulink的实时仿真能力与MATLAB的分析功能形成完美互补。

构建一个电机速度控制系统的Simulink模型（motor_control.slx），包含以下关键模块：

• 电机传递函数：1/(Js + b)，其中J=0.01 kg·m²，b=0.1 N·m·s
• PID控制器：比例增益Kp=2，积分时间Ti=0.5s
• 噪声注入：带宽受限白噪声模拟传感器干扰

频率响应测试配置步骤：

1. 在Simulink模型中插入Frequency Response分析模块
2. 设置扫频范围：0.1Hz到100Hz，对数间隔
3. 指定输入输出端口（激励信号与测量点）
4. 运行频率扫描并自动生成Bode图

% 从Simulink模型导出频率响应数据 load_system('motor_control.slx') in = 'motor_control/Reference'; out = 'motor_control/Speed'; op = operpoint('motor_control'); % 创建频率响应对象 frd_model = frestimate('motor_control', op, in, out, logspace(-1,2,200)); % 对比理论模型与仿真结果 figure bode(frd_model, 'r', G_theoretical, 'b--') legend('Simulink实测', '理论模型') title('模型验证对比')

Simulink频率分析优势：

• 非线性系统处理：当系统包含饱和、死区等非线性环节时仍能有效分析
• 噪声影响评估：可量化测量噪声对频率特性的影响程度
• 控制器验证：在实际闭环条件下测试控制器性能
• 硬件在环(HIL)：支持连接实际控制器硬件进行测试

对于数字控制系统，Simulink还能自动处理离散化效应。设置采样时间为0.01s后，系统会自动考虑零阶保持器(ZOH)引入的相位滞后：

% 离散系统频率分析 G_digital = c2d(G, 0.01, 'zoh'); figure bode(G, G_digital) legend('连续系统', '离散系统(Ts=0.01s)') title('采样对频率特性的影响')

这种分析对数字控制器设计至关重要，它能准确预测采样引起的相位损失，避免实际系统中出现意外的稳定性问题。

5. 工程案例分析：伺服系统稳定性优化实战

某工业机械臂关节伺服系统出现低频振荡问题，原始设计开环传递函数为： $$ G(s) = \frac{200}{s(0.5s+1)(0.02s+1)} $$

通过MATLAB分析发现相位裕度不足（仅35°），需进行控制器优化。采用相位超前补偿： $$ C(s) = \frac{0.1s + 1}{0.025s + 1} $$

优化步骤实现：

% 原始系统分析 G = tf(200, conv([1 0], conv([0.5 1], [0.02 1]))); margin(G) % 显示原始稳定裕度 % 设计补偿器 zc = 1; % 零点频率1 rad/s pc = 40; % 极点频率40 rad/s C = tf([1/zc 1], [1/pc 1]); % 补偿后系统分析 figure margin(C*G) % 显示补偿后稳定裕度 % 时域验证 t = 0:0.001:2; figure step(feedback(G, 1), 'b') % 原始系统阶跃响应 hold on step(feedback(C*G, 1), 'r') % 补偿后响应 legend('原始系统', '补偿后系统') title('时域性能对比')

优化效果对比表：

这个案例展示了频率法在工程调试中的实际价值——通过Bode图分析定位问题根源，设计针对性的补偿网络，最终实现系统性能的全面提升。整个过程将理论分析、软件工具和工程经验紧密结合，体现了现代控制工程师的典型工作流程。

6. 高级技巧：自动化报告生成与批量处理

对于需要分析大量系统变体的工程团队，手动处理每个案例效率低下。MATLAB支持脚本化批量处理和自动报告生成，极大提升工作效率。

批量分析示例：

% 定义多个系统配置 systems = { tf(100, [1 5 100]), % 二阶系统 tf(50, conv([1 0], [1 10])), % 含积分环节 tf([1 1], conv([1 0 0], [0.1 1])) % 非最小相位 }; % 创建分析报告 report = ['<html><body><h1>系统频率特性分析报告</h1>', date]; for i = 1:length(systems) G = systems{i}; % 生成分析图形 fig = figure('Visible', 'off'); subplot(2,1,1) bode(G) title(['系统', num2str(i), ' Bode图']) subplot(2,1,2) nyquist(G) title(['系统', num2str(i), ' Nyquist图']) % 保存图形到报告 imgfile = ['system', num2str(i), '.png']; saveas(fig, imgfile) report = [report, '<h2>系统', num2str(i), '</h2>', ... '<img src="', imgfile, '" width="800">']; % 计算性能指标 [Gm,Pm,Wcg,Wcp] = margin(G); report = [report, '<ul>', ... '<li>增益裕度: ', num2str(20*log10(Gm)), ' dB</li>', ... '<li>相位裕度: ', num2str(Pm), ' deg</li>', ... '<li>穿越频率: ', num2str(Wcp), ' rad/s</li></ul>']; end % 保存报告 report = [report, '</body></html>']; fid = fopen('analysis_report.html', 'w'); fprintf(fid, report); fclose(fid);

这种自动化流程特别适合产品系列化开发或参数敏感性研究。工程师可以一次性提交多个设计方案，通过脚本自动生成包含所有关键指标的比较报告，大幅缩短设计迭代周期。

对于更复杂的工程系统，还可以结合MATLAB的面向对象编程特性构建专用分析类：

classdef FrequencyAnalyzer properties System MarginInfo FrequencyResponse end methods function obj = FrequencyAnalyzer(sys) obj.System = sys; [Gm,Pm,Wcg,Wcp] = margin(sys); obj.MarginInfo = struct('GainMargin',Gm, ... 'PhaseMargin',Pm, ... 'CrossoverFreq',Wcp); obj.FrequencyResponse = freqresp(sys, logspace(-2,3,500)); end function plotComparativeBode(obj, sys2) figure bode(obj.System, 'b', sys2, 'r--') legend('原始系统', '对比系统') end end end % 使用示例 G1 = tf(1, [1 1 0]); analyzer = FrequencyAnalyzer(G1); G2 = tf(2, [1 1 0]); analyzer.plotComparativeBode(G2)

这种面向对象的设计模式使分析工具更具扩展性和复用性，特别适合构建企业内部的专用分析平台。