简介：面向无线通信、信号处理方向的学习者与研究人员，这份MATLAB代码资源聚焦多输入多输出MIMO系统的广义奈奎斯特曲线绘制，帮助评估系统稳定性与稳定边界。资源包共7个文件，全部为.m脚本，压缩包约5KB，涵盖传递函数符号化转换、模型降阶、广义奈奎斯特绘图及示例主程序等模块，结构紧凑、便于逐文件阅读与二次修改。已有3163人学习下载，说明其在教学与科研场景中具备一定参考价值。读者可借助示例脚本直接运行，观察MIMO系统在复平面上的稳定边界，并结合符号运算与降阶工具理解多输入多输出耦合下的频率响应与相位特性，从而把通信理论与MATLAB编程实践结合起来，适用于课程实验、论文复现与系统分析入门。

1. 多输入多输出系统的稳定性判据，为什么单回路奈奎斯特不够用

做 MIMO 控制系统频域分析的人，迟早会撞上一个尴尬局面：每个回路单独画奈奎斯特曲线都稳稳当当，闭环跑起来却发散。这不是玄学，是单回路判据在多输入多输出系统里失效了。单回路奈奎斯特只回答“某个回路增益裕度够不够”，它看不见回路之间的耦合——一路的相位滞后会通过交叉通道把另一路推向不稳定。广义奈奎斯特（Generalized Nyquist）判据正是为这件事准备的：把开环传递矩阵 L(s)=G(s)K(s) 的特征值轨迹画在复平面上，用特征轨迹包围临界点的圈数判断闭环稳定性。这份 MATLAB 资源做的就是这件事——给定 MIMO 系统的传递矩阵或状态空间模型，自动扫频、求特征值、绘制广义奈奎斯特曲线并给出稳定性结论。适合正在做多变量频域设计、需要验证解耦控制器或鲁棒控制器闭环稳定性的从业者，也适合课程设计阶段需要把判据跑通、把图画出规范样子的同学。

2. 广义奈奎斯特判据的数学骨架：从传递矩阵到特征轨迹

2.1 判据本身：特征值轨迹与临界点的包围关系

先把判据说清楚，否则代码只是黑匣子。对单位负反馈闭环，开环传递矩阵为 L(s)，闭环特征多项式与 det(I+L(s)) 直接相关。广义奈奎斯特判据的表述是：令 s 沿奈奎斯特围线（虚轴加右半平面大圆弧）走一圈，L(s) 的特征值 λ_i(s) 在复平面上画出若干条轨迹，这些轨迹对临界点 (-1, 0) 的净包围圈数 N，满足 Z = N + P，其中 P 是 L(s) 在右半平面的极点数，Z 是闭环右半平面极点数。Z=0 即闭环稳定。

这里有两个容易翻车的细节。第一，包围圈数要按“净”算，逆时针为正、顺时针为负，方向搞反结论直接颠倒。第二，P 必须从开环模型里数出来，不能默认是 0——很多实际系统开环就带右半平面极点（比如某些欠阻尼对象或不稳定平台），P 不为零时判据的结论会完全不同。资源里如果只画图不给 P 和 N 的计数，那图只能看个热闹。

单回路判据是这套东西在 L(s) 退化为标量时的特例，此时特征值就是 L(s) 本身，轨迹只有一条。MIMO 下特征值有多条，每条都要数包围，且特征值轨迹在特征值重合点附近会出现“交换”现象，排序错乱会让曲线看起来断裂——这是绘图层面最常见的视觉坑，后面避坑章节会专门讲。

2.2 为什么用特征值而不是行列式

理论上 det(I+L(s)) 的幅角变化也能判稳，但工程上没人这么画。原因很实际：行列式把多个回路的耦合压成一个复数，曲线一旦穿过临界区域，你根本不知道是哪个通道、哪个频段出的问题；而特征值轨迹保留了通道信息，哪条轨迹贴近 (-1,0)、在什么频率贴近，一眼能看出来，方便回头调增益或加补偿。代价是特征值分解在数值上对重特征值敏感，扫频点稀疏时轨迹会跳变，所以扫频密度和特征值排序策略是这份资源能不能用的关键。

我一般会要求扫频点在关键频段（穿越频率附近）加密，至少每十倍频程 200 点起步，特征值排序用“相邻频率点最近邻匹配”而不是 MATLAB 默认的 eig 输出顺序。默认顺序在特征值交叉处会随机交换，画出来的曲线像被剪断的绳子。

2.3 从模型到特征轨迹：核心计算流程

不管资源封装成函数还是脚本，底层流程是固定的四步：构造频率向量、逐点计算 L(jω)、逐点做特征值分解、按最近邻原则重排特征值序列后绘图并计数。下面这段是我按这个流程写的可复现骨架，参数含义在代码后说明。

% 广义奈奎斯特曲线绘制骨架 % G: 开环传递矩阵（tf 或 ss 对象），n×n % w: 扫频向量，单位 rad/s function [lam, w] = gen_nyquist_traj(G, w) n = size(G, 1); lam = zeros(n, numel(w)); % 每行存一条特征值轨迹 prev = []; % 上一频点的特征值，用于最近邻排序 for k = 1:numel(w) Lk = evalfr(G, 1j*w(k)); % 计算 L(jw) e = eig(Lk); % 特征值分解 if isempty(prev) e = sort(e); % 首点按实部排序定基准 else % 最近邻匹配：每个新特征值找上一频点最近的那个 idx = zeros(n,1); used = false(n,1); for i = 1:n d = abs(e - prev(i)); d(used) = inf; [~, idx(i)] = min(d); used(idx(i)) = true; end e = e(idx); end lam(:,k) = e.'; prev = e; end end

逻辑说明：evalfr逐点求值避免 tf 对象在向量化调用时的维度陷阱；eig得到该频点的 n 个特征值；最近邻匹配保证第 i 条轨迹在相邻频点间连续，这是曲线不断裂的前提。参数上，w建议用logspace(log10(wmin), log10(wmax), N)生成，wmin 取最低关注频率的 1/10，wmax 取最高关注频率的 10 倍，N 在 2000 以上；G若是状态空间模型，evalfr同样适用，但要注意ss模型的 D 矩阵非零时高频段特征值会趋于 D 的特征值，绘图范围要留够。

画图时把lam的每一行plot(real(lam(i,:)), imag(lam(i,:)))，再叠加plot(-1, 0, 'r+')标出临界点，最后用axis equal保证包围圈数在视觉上不被拉伸误导——这一条很多人忽略，坐标轴不等比时“看起来包围了”和“实际包围了”是两回事。

3. 在 MATLAB 里把曲线跑出来：模型构造、扫频与稳定性计数

3.1 构造一个可验证的 MIMO 开环模型

拿一个 2×2 的传递矩阵做例子，通道间带耦合，这样能体现广义奈奎斯特和单回路的差别。下面用tf直接搭，交叉项故意给一个相位滞后的耦合。

s = tf('s'); % 2x2 开环传递矩阵，对角为主通道，非对角为耦合通道 G11 = 10/(s*(s+1)); G12 = 0.8/(s+2); G21 = 0.5/(s+3); G22 = 6/((s+0.5)*(s+2)); G = [G11 G12; G21 G22]; % 控制器（示例：对角比例+积分） K = diag([1.2*(1+1/(0.5*s)), 0.9*(1+1/(0.8*s))]); L = G*K; % 开环传递矩阵 P = count_rhp_poles(L); % 数右半平面极点，见 3.3

逻辑说明：G用tf拼成矩阵，MATLAB 会自动处理 MIMO 传递矩阵的运算；K用diag构造对角控制器，这是最常见的解耦前配置。参数上，主通道增益和积分时间决定穿越频率，耦合项系数 0.8 和 0.5 控制通道间影响强度——把这两个数调大，你会看到特征轨迹明显向 (-1,0) 靠拢，这就是耦合把系统推向不稳定的直观表现。count_rhp_poles不是内置函数，需要自己写或用pole加实部判断，见下一节。

3.2 扫频与特征轨迹计算

把 2.3 的骨架接上，扫频范围按系统带宽定。这个例子里主通道穿越频率在 3 rad/s 附近，扫频取 0.01 到 1000 rad/s。

w = logspace(-2, 3, 3000); % 0.01 ~ 1000 rad/s，3000 点 [lam, w] = gen_nyquist_traj(L, w); % 复用 2.3 的函数 figure; hold on; grid on; for i = 1:size(lam,1) plot(real(lam(i,:)), imag(lam(i,:)), 'LineWidth', 1.2); end plot(-1, 0, 'r+', 'MarkerSize', 12, 'LineWidth', 2); axis equal; % 关键：等比例坐标 xlabel('Re'); ylabel('Im'); title('Generalized Nyquist: eigenvalue loci of L(jw)');

逻辑说明：logspace保证低频和高频都有足够分辨率；循环逐条画特征轨迹；红色十字标临界点；axis equal是判读包围关系的前提。参数上，3000 点是经验值，若曲线在穿越区仍有折角，加到 5000；w的下限不要低于系统最低极点实部的 1/100，否则低频段轨迹挤成一团看不出包围方向。

3.3 数 P 和 N：把结论从图里算出来

图能看，但结论要算。P 用pole求开环极点后数实部为正的个数；N 用轨迹对 (-1,0) 的环绕数，工程上常用“轨迹穿越负实轴 (-∞,-1) 区间的次数与方向”来累计，或者直接用幅角增量法。

function P = count_rhp_poles(L) p = pole(L); % 所有极点（含重根展开） P = sum(real(p) > 1e-8); % 实部为正计为右半平面 end function N = winding_number(lam, w, crit) % 幅角增量法：对每条轨迹累加相对临界点的角度变化 N = 0; for i = 1:size(lam,1) z = lam(i,:) - crit; dphi = diff(unwrap(angle(z))); % 展开相位避免跳变 N = N + sum(dphi)/(2*pi); % 净圈数，逆时针为正 end end

逻辑说明：count_rhp_poles用pole拿全部极点，实部阈值 1e-8 是为了滤掉数值零；winding_number用unwrap展开相位后累加角度增量，除以 2π 得净圈数。参数上，crit传-1（即 -1+0j）；unwrap对扫频点密度敏感，点太稀时相位跳变超过 π 会误判，这也是为什么前面强调 3000 点起步。最后Z = N + P，Z 为 0 判稳。注意这里的 N 是净圈数，正负号取决于扫频方向，扫频从低频到高频时逆时针为正，若你的实现扫频反向，符号要取反再代入。

3.4 和单回路结果对照，确认判据真的在起作用

把L的对角元单独拿出来各画一条单回路奈奎斯特，和特征轨迹放一起看。这个例子里如果耦合系数调到 1.5 以上，单回路两条曲线都离 (-1,0) 有距离，但特征轨迹已经绕进临界点——这就是开头说的“单回路都稳、闭环发散”。做这一步对照不是为了好看，是为了在交付前确认你的广义奈奎斯特实现没有把耦合项漏掉。常见错误是误用G(:,:,1)之类的索引把 MIMO 模型降成单入单出，图是画出来了，判据却退化了。

4. 避坑与排查：特征轨迹绘制的五类翻车现场

4.1 曲线在特征值交叉处断裂

现象：图上某条轨迹突然跳到另一条的位置，看起来像断成两截。原因：eig返回的特征值顺序在每个频点不保证连续，交叉点附近顺序随机交换。解决：用 2.3 的最近邻匹配重排，且首点用sort定基准；若仍有跳变，把该频段扫频点加密到 5000 以上再试。

4.2 包围圈数算出来和肉眼不一致

现象：图上明明绕了临界点，winding_number返回 0 或符号相反。原因：一是坐标轴没axis equal，视觉包围是拉伸假象；二是unwrap在稀疏频点处相位跳变超过 π 被误展开；三是扫频方向与符号约定不一致。解决：先axis equal复核视觉，再把扫频点加密，最后核对扫频是从低频到高频还是反向，反向时 N 取负。

4.3 P 默认按 0 处理导致结论错误

现象：开环明明有右半平面极点，判据却按 P=0 算，得出“稳定”的错误结论。原因：直接套 Z=N 而没数开环极点，或pole对 MIMO 模型返回的是所有通道极点拼接，重根被重复计数。解决：用count_rhp_poles显式数，重根用uniquetol去重后再数；对ss模型确认pole返回的是系统极点而非单通道极点。

4.4 高频段轨迹发散到画布外

现象：曲线在高频处冲到复平面很远，图的比例被拉坏，临界点附近挤成一团看不清。原因：ss模型 D 矩阵非零或传递函数分子阶次不低于分母时，高频特征值趋于常数或发散，扫频上限设太高。解决：扫频上限收到最高关注频率的 10 倍以内，绘图时用xlim/ylim裁到临界点附近区域，或对轨迹做分段绘制。

4.5 把 MIMO 模型误降成单回路

现象：特征轨迹只有一条，或条数与系统维数不符。原因：模型构造时用了G(1,1)之类的单元素索引，或evalfr传入了标量化的模型。解决：绘图前size(L)确认是 n×n，evalfr返回矩阵后size再确认一次；特征值条数应等于系统输入维数。

5. 进阶：把广义奈奎斯特接进鲁棒性验证与批量扫参

单次画图只是起点，真正省时间的是把它做成可批量跑的验证函数。我一般会包一层扫参循环，对控制器增益或耦合系数做网格扫描，每个点算 Z，输出稳定域图。下面这段是骨架，gain_grid是待扫的增益向量，Zmap存每个组合的闭环右半平面极点数。

gains = 0.5:0.1:2.5; Zmap = zeros(numel(gains), numel(gains)); w = logspace(-2, 3, 2000); for a = 1:numel(gains) for b = 1:numel(gains) K = diag([gains(a)*(1+1/(0.5*s)), gains(b)*(1+1/(0.8*s))]); L = G*K; [lam, ~] = gen_nyquist_traj(L, w); P = count_rhp_poles(L); N = winding_number(lam, w, -1); Zmap(a,b) = N + P; % 0 为稳定 end end imagesc(gains, gains, Zmap == 0); % 稳定域可视化

逻辑说明：双层循环扫两个通道的增益，每个组合重算特征轨迹、P、N，Zmap==0标出稳定区。参数上，扫参步长 0.1 是精度和耗时的折中，2000 点扫频在 25×25 网格下大约几十秒量级，若嫌慢可先粗扫再在边界附近细扫。这张稳定域图比单张奈奎斯特曲线有用得多——它直接告诉你两个通道增益的可行组合边界，而单回路裕度给不出这个信息。

验证方法上，我习惯用闭环极点做交叉校验：对稳定域内的点，用feedback(L, eye(2))求闭环极点，确认全部实部为负；对边界点，闭环极点应恰好贴近虚轴。这一步是后悔药——广义奈奎斯特的计数一旦有符号或密度问题，闭环极点会立刻暴露，比盯着曲线数圈数可靠。另外，特征值轨迹在重特征值附近对参数极敏感，扫参时若发现某网格点 Z 突变，先怀疑是数值问题而不是系统真不稳定，把该点扫频加密重算一遍再下结论。

从那以后我每次交付频域稳定性结论，都强制走一遍“特征轨迹图 + P/N 显式计数 + 闭环极点交叉校验”三件套，少一步都不敢签字。希望帮到你。

本文还有配套的精品资源，点击获取