简介：一套面向无线通信与自动控制专业学生及研究者的MATLAB工具包，围绕多输入多输出（MIMO）系统广义奈奎斯特（Nyquist）曲线绘制而设计，可用于分析多变量系统的闭环稳定性与奈奎斯特判据在MIMO场景下的扩展应用。资源以RAR压缩包形式提供，共7个文件且均为.m脚本，压缩包仅5KB，内容涵盖系统状态空间转换、符号函数处理、模型降阶以及广义奈奎斯特曲线的主绘制函数等，覆盖从系统建模到曲线绘制的完整流程。已有3162人浏览学习。借助脚本示例与主绘制函数，读者可直接运行观察不同传递函数矩阵下广义奈奎斯特曲线的形态，结合理论分析各输入输出通道间的相位与幅值特性，快速掌握此类稳定性分析工具的使用方法。适合正在学习MIMO系统课程、开展无线通信稳定性研究或需要完成相关实验报告的读者参考使用。
1. 多输入多输出系统的广义奈奎斯特曲线：从标量到矩阵的一次跨越
做过多变量控制系统的人基本都吃过亏：明明每个单回路都稳定，联调起来却莫名其妙振荡甚至发散。原因在于传统奈奎斯特判据只适用于单输入单输出（SISO）系统，而MIMO系统的回比矩阵是复数矩阵，行列式的零点无法直接用标量奈奎斯特图判断。这套MATLAB资源解决的正是这个问题——它帮你用广义奈奎斯特曲线（特征轨迹法）绘制MIMO系统的回比矩阵特征值轨迹，把矩阵判据拆成一组标量曲线来分析，从而判断闭环系统的稳定性。适合做多变量鲁棒控制、解耦控制或者正在为课程设计发愁的从业者，尤其是被多变量频域分析卡住、想直接拿到可运行脚本的那些人。接下来我会从数学原理讲起，再给出完整的MATLAB实现方案、参数设置和几个容易翻车的细节。
2. 广义奈奎斯特判据：为什么特征轨迹能替代行列式判据
2.1 从回比矩阵到特征值映射的推演过程
MIMO系统的稳定性分析绕不开回比矩阵L(s)，它的定义是开环传递函数矩阵G(s)与控制器矩阵K(s)的乘积，即L(s) = G(s)K(s)。对于SISO系统，我们画出L(jω)的奈奎斯特曲线，看它是否包围临界点(-1, j0)就能判断闭环稳定性。但MIMO系统的L(jω)是一个m×m的复数矩阵，直接在同一个复平面上画矩阵元素是没意义的，因为矩阵的行列式零点等价于特征值穿越原点的次数，但特征值本身是复数函数。
广义奈奎斯特判据的核心思想是把矩阵的行列式det(I + L(s))的零点问题转化为L(s)特征值轨迹的包围问题。对于每个频率ω，L(jω)都有m个特征值λi(jω)，其中i = 1, 2, ..., m。当这些特征值轨迹在复平面上绕(-1, j0)点的净包围次数等于开环右半平面极点数时，闭环系统稳定。这个结论在文献中叫广义奈奎斯特判据，它是经典奈奎斯特判据在矩阵情形的自然推广。
需要特别注意的是，特征值轨迹必须逐根追踪，不能把不同频率下的特征值打乱顺序。L(jω)的特征值在数值上可能发生交叉，比如某些频率下第1个特征值和第2个特征值的模很接近，数值求解器可能会交换它们的顺序，导致轨迹看起来有突跳。解决方法是利用特征值的光滑性进行排序或连续性追踪，后面实现部分我会给出具体方案。
2.2 奈奎斯特阵列与Gershgorin带的关系解读
广义奈奎斯特曲线还有一种可视化补充方式叫奈奎斯特阵列。它的做法是把L(jω)的对角元素作为主曲线，然后用Gershgorin圆盘定理给每条对角曲线加一个半径带。Gershgorin定理说矩阵的每个特征值都位于以对角元素为中心、半径为该行（或列）非对角元素绝对值之和的圆盘内。这个性质用来做对角优势判断很有用——如果每个Gershgorin带都不包含(-1, j0)点，就能直接得出闭环稳定性的结论，而无需逐点求解特征值。
这个资源的实现里同时绘制了特征轨迹和Gershgorin带，两者互为验证。特征轨迹是精确结果，Gershgorin带是保守估计。实际使用中，如果某个Gershgorin带包住了临界点但特征轨迹没有，说明系统的对角优势不够强，这时候即使闭环稳定，鲁棒性也可能很差。我在做某跨平台系统的鲁棒性分析时遇到过类似情况，后来改成先画Gershgorin带粗判断，再用特征轨迹精确判断，分析效率提高不少。
2.3 特征轨迹法的适用边界与前提条件
特征轨迹法不是万能的，它要求开环传递函数矩阵L(s)在虚轴上没有极点，如果虚轴上有极点或非常接近虚轴的极点，直接代入s = jω计算会得到异常大的数值，导致曲线飞出坐标范围。通常的预处理是先把虚轴上的极点稍微向左平移（极点配置到s = -ε，比如ε取1e-3），或者在频点序列中跳过这些极点的谐振峰。
另一个前提是L(s)必须是方阵，也就是输入输出维数相等。对于非方阵MIMO系统，比如2输入3输出，回比矩阵是3×3方阵（输出维数），但控制器的参数化会更复杂。这个资源的脚本默认处理方阵情形，非方阵需要自行扩展特征值分解的维度，我后面会补充一个处理思路。
3. MATLAB实现广义奈奎斯特曲线：核心函数与完整流程
3.1 频响矩阵的批量计算方法
绘制广义奈奎斯特曲线的第一步是计算L(jω)在预设频点上的复数矩阵值。直接写嵌套循环逐个频点求解确实简单，但MIMO系统状态空间维度高时速度极慢。这个资源默认采用状态空间模型批量计算频响，避免了频繁调用函数带来的开销，具体实现代码如下。
function Lfreq = compute_loop_freq(A, B, C, D, K, freqs) % 批量计算回比矩阵L(jw) = C*(jwI-A)^(-1)*B*K + D*K % 输入: A,B,C,D - 被控对象状态空间矩阵 % K - 控制器增益矩阵 (m x m) % freqs - 频率点向量 (rad/s)，建议对数均匀分布 % 输出: Lfreq - 三维数组, Lfreq(:,:,k) 表示第k个频点的L(jw) n_freq = length(freqs); [~, m] = size(B); % B的列数 = 输入通道数 Lfreq = zeros(m, m, n_freq); I = eye(size(A)); for k = 1:n_freq w = freqs(k); sI_A = 1i * w * I - A; % 用左除求解频响矩阵，避免显式求逆 Gjw = C * (sI_A \ B) + D; Lfreq(:,:,k) = Gjw * K; end end这段代码的核心思路是先用左除对(sI - A)的线性系统进行求解，而不是直接对矩阵求逆，这样数值稳定性更好。参数上最需要注意的是A矩阵可能很大，直接对高维矩阵用左除会消耗不少内存，所以我一般会在调用前把状态空间模型先做一次最小实现（minreal）处理，把不可控不可观的维度削掉再传入。另外，计算频响时用的是1i * w，复数单位不能写成i，否则会被当成变量，这是MATLAB新人经常踩的坑。
3.2 特征值轨迹逐根排序与绘制
得到三维数组Lfreq之后，核心工作是在每个频点对矩阵做特征值分解，然后把特征值按连续性排序并绘制。这里的难点是特征值排序，如果直接调用eig函数，相邻频点返回的特征值顺序可能交换，绘制出来的轨迹会有锯齿状跳变，看起来就像两条曲线互相穿越。下面的代码用最简单的最小距离匹配策略做排序。
function plot_generalized_nyquist(Lfreq, freqs) % 绘制广义奈奎斯特曲线的特征轨迹 % 策略: 逐频点求解特征值, 按与前一个频点的欧氏距离最小进行匹配排序 n_freq = size(Lfreq, 3); m = size(Lfreq, 1); colors = lines(m); % 为每条轨迹分配颜色 figure; hold on; grid on; axis equal; % 保持实部虚部比例一致 prev_evals = []; for k = 1:n_freq evals = eig(Lfreq(:,:,k)); % 当前频点的m个特征值 if k > 1 % 匹配排序: 找使总距离最小的排列(简化为贪婪策略) for idx = 1:m dists = abs(evals - prev_evals(idx)); [~, match] = min(dists); % 交换到对应位置(此处为示意, 实际建议用匈牙利算法) tmp = evals(idx); evals(idx) = evals(match); evals(match) = tmp; end end prev_evals = evals; for idx = 1:m plot(real(evals(idx)), imag(evals(idx)), ... 'o', 'Color', colors(idx,:), 'MarkerSize', 4); if k > 1 % 连接前后点, 形成连续轨迹 plot([real(prev_evals(idx)) real(evals(idx))], ... [imag(prev_evals(idx)) imag(evals(idx))], ... '-', 'Color', colors(idx,:), 'LineWidth', 1.2); end end prev_evals = evals; end plot(-1, 0, 'kx', 'MarkerSize', 12, 'LineWidth', 2); xlabel('Re'); ylabel('Im'); title('Generalized Nyquist Curves (Characteristic Loci)'); end这段排序代码用了贪心策略，实现简单但并非最优匹配。实际工程中频点足够密时，贪心策略基本够用；频点稀疏或特征值轨迹交叉频繁时，建议用匈牙利算法做全局最优匹配。另一个隐藏问题是首频点特征值的初始顺序是任意的，所以每条轨迹的颜色在起点处可能和理论上的第几通道不严格对应，这并不影响稳定性判断，但如果需要把轨迹对应到具体输入输出通道，就要在首频点做一次人工校准。
3.3 完整脚本：从状态空间模型到输出曲线的串联流程
资源包里的主脚本把前面两个函数串联起来，同时附加了Gershgorin带的绘制功能。它从给定的状态空间模型出发，自动选择频率范围，计算频响，绘制特征轨迹和Gershgorin带，并把临界点(-1, j0)标出来。主脚本的流程如下。
% main_mimo_nyquist.m % 广义奈奎斯特曲线绘制主脚本 % 示例: 2x2 MIMO系统, 状态空间模型由文件demo_model.m生成 clear; clc; close all; % 步骤1: 加载或定义状态空间模型 (2输入2输出, 4阶) % 这里直接用随机矩阵构造一个稳定的开环对象示范 rng(2024); A = [-2 0.5 0 0; -0.4 -1.5 0.8 0; 0 -0.6 -3 0.2; 0.5 0 -0.3 -2.8]; B = [1 0; 0 1; 0.5 0.3; 0 0.6]; C = [1 0.2 0 0.4; 0 0.3 1 0]; D = zeros(2,2); % 被控对象无直通项 K = [2.5 0.8; -0.5 2.0]; % 控制器增益矩阵 % 步骤2: 生成对数均匀分布的频点 freqs = logspace(-2, 2, 800); % 0.01 rad/s 到 100 rad/s % 步骤3: 批量计算回比矩阵频响 Lfreq = compute_loop_freq(A, B, C, D, K, freqs); % 步骤4: 绘制广义奈奎斯特特征轨迹 plot_generalized_nyquist(Lfreq, freqs); % 步骤5: 计算并叠加Gershgorin带(简化实现) figure; hold on; grid on; axis equal; for k = 1:length(freqs) Lk = Lfreq(:,:,k); r = sum(abs(Lk), 2) - abs(diag(Lk)); % 行Gershgorin半径 for idx = 1:size(Lk,1) center = Lk(idx,idx); theta = linspace(0, 2*pi, 60); plot(real(center) + r(idx)*cos(theta), imag(center) + r(idx)*sin(theta), '-'); end end plot(-1, 0, 'rx', 'MarkerSize', 14);主脚本里的参数矩阵是随机生成的，实际使用时要替换成你自己的被控对象模型和控制器参数。频率范围选择上，logspace(-2, 2, 800)表示从0.01到100 rad/s取800个点，这个密度对大多数工程问题是足够的。如果系统有更低的谐振频率，可以把下限向下扩展，但要注意频点增多会显著增加特征值分解的计算量。
3.4 自定义频响模型：直接输入传递函数矩阵的替代方案
很多场景下你没有状态空间模型，只有传递函数矩阵G(s)的数值数据，或者直接从频域辨识得到频响数据。这种情况下没必要强行构造状态空间，可以直接把频响数据组织成三维数组传入绘制函数。下面给出一个替代入口，它接收的是频响数据而非状态空间矩阵。
function plot_nyquist_from_frd(frd_data, K, freqs) % 从频响数据直接绘制广义奈奎斯特曲线 % frd_data: m x m x n_freq 复数数组, 表示不同频率下的G(jw) % K : 控制器增益矩阵 n_freq = size(frd_data, 3); m = size(frd_data, 1); Lfreq = zeros(m, m, n_freq); for k = 1:n_freq Lfreq(:,:,k) = frd_data(:,:,k) * K; end plot_generalized_nyquist(Lfreq, freqs); end这种替代方案非常适合控制系统工具箱里有现成频响数据但不想做系统辨识的场景。它还避免了状态空间模型维数过高导致左除求解慢的问题，直接把频响数据喂进来就行。需要注意频率向量freqs必须和frd_data的第三维一一对应，长度不一致时后面的曲线会整体错位，我见过有人把单位搞混（Hz和rad/s混用）导致曲线形态完全对不上的案例。
4. 避坑与常见问题排查：特征轨迹绘制中的五个典型翻车点
4.1 特征值交叉与轨迹跳变问题
现象：绘制的特征轨迹在某段频率区间出现折线拐角，或者两条曲线像被强行打断了，形成毛刺状连接。
原因：L(jω)的两个特征值在某个频率附近数值上非常接近，eig函数的返回顺序随矩阵数值扰动发生交换，导致前后频点的同一条轨迹断开了。
解决：把频点加密到原来的3到5倍是最直接的办法。频点间隔变密后，特征值变化的位移量变小，贪心匹配的失败率会显著下降。如果加密后仍然出现跳变，就要用匈牙利算法做全局匹配，MATLAB中可以调用matchpairs函数。脚本里默认的800个频点对大多数系统已经够用，但系统阻尼比极低、谐振峰尖锐时需要临时加密。
4.2 虚轴极点导致曲线飞出坐标范围
现象：曲线在某个频率附近突然发射出很长的线段，直达坐标框外，整个图的比例尺被拉伸，其余特征轨迹被压成一条细线。
原因：开环传递函数在虚轴上存在极点。代入s = jω时，如果ω正好等于极点的虚部，频响矩阵的某个元素趋向无穷大，特征值也随之发散。
解决：先对开环系统做极点分析，找出实部为零或实部绝对值小于1e-6的极点，然后把系统稍微稳定化处理——把极点向左平移ε即可，常见的做法是给积分器加一个小的阻尼项。平移后系统的低频特性有细微变化，但对稳定性判断的影响在工程允许范围内。另一个办法是人为跳过谐振峰附近的频点，但这个做法会掩盖真实的穿越信息，不是推荐方案。
4.3 Gershgorin带计算时误用列半径与行半径
现象：绘制的Gershgorin带宽度和特征轨迹的覆盖范围对不上，看起来带子特别宽或特别窄，最终判据结论不可信。
原因：Gershgorin定理有两种形式，行形式和列形式，它们分别给出特征值的不同的包含区域。如果混用了行半径计算对角元素中心、列半径计算宽度，得到的带子就不是有效的包含区域。
解决：统一使用行半径，即r_i = Σ_{j≠i} |L_{ij}|，中心取L_{ii}。实际判断时我常把行半径和列半径各画一组，取两者的交集作为特征值的实际包含区域，能明显收紧保守性。不过要注意行、列两组带子分别用不同的线型或颜色区分，否则图面会很凌乱。
4.4 临界点包围次数判读时忽略了开环不稳定极点
现象：画出曲线后，凭直觉数了绕(-1, j0)点的圈数，然后直接下结论闭环稳定，结果和实际仿真不符。
原因：广义奈奎斯特判据要求包围次数等于开环右半平面极点数P，稳定条件是N = P。很多人数出包围次数却不检查P，或者把逆时针包围当作正方向，导致正负号弄反。
解决：先对开环系统L(s)计算右半平面极点数，在MATLAB中用pole(G*K)然后检查实部大于0的个数。然后约定逆时针包围为正。最稳妥的验证方式是画完广义奈奎斯特曲线后，再用闭环特征值直接计算一次闭环极点，如果极点实部全为负而判据结果说稳定，两者互相印证即可。
4.5 复数单位i与循环变量i的命名冲突
现象：脚本报错提示找不到变量i或变量被当作复数单位使用，频响计算结果全部为零或报维度错误。
原因：MATLAB中i和j默认是复数单位，但很多人习惯用i做循环索引，在循环中执行了i = 1这样的赋值后，后续代码再写1ijw时i被重新赋值，导致复数计算全部出错。
解决：养成两个习惯，一是循环变量尽量用k、idx这类字母，把i和j保留给复数单位；二是在频响计算代码中统一使用1i而不是i。脚本里我全部用了k做循环变量就是为了避开这个坑。如果你接手别人的代码，先全局搜索是否存在对i或j的赋值语句，存在的话优先改掉。
5. 参数调优与仿真效率：频点选择与算法加速手段
5.1 频点分布策略：线性、对数与自适应加密
绘制广义奈奎斯特曲线的频点分布直接决定曲线的平滑度和计算量。线性分布适合带宽窄、谐振峰平坦的系统，但MIMO系统的动态范围往往很宽，低频段和高频段的增益差异可能达到几个数量级，用线性分布会导致低频段的曲线点过于稀疏，看不清穿越细节。
对数分布是多数情况下的首选，logspace生成的点在低频频段密、高频频段疏，符合工程系统频响的分布规律。但对数分布也有盲区，如果系统在某个较高频段存在尖锐谐振峰，对数分布的频点间隔可能跨过峰顶，导致谐振尖峰被采样模糊掉。面对这种情况，我一般在已知谐振频率附近做一次线性加密，把logspace和linspace两段拼接到一起，效果比单纯加密整个频带更好。
自适应加密是更高级的做法：第一次粗扫出特征值轨迹后，定位相邻两点之间梯度变化最大的频段，在这些频段插入新频点，迭代两到三轮。资源里提供了一个参数n_refine控制加密轮数，默认设为2，对于典型的多变量飞行器模型，加密后曲线的平滑度和稳定性判定精度都能明显提升。
5.2 高维MIMO系统的计算瓶颈与降维预处理
当系统维数m上升到5阶以上时，逐频点特征值分解的耗时开始变得不可忽略。考虑一个10输入10输出的系统，800个频点就要做800次10×10矩阵的特征值分解，虽然单次耗时只有毫秒级，但加上频响矩阵的求解后总耗时可能达到数十秒。更关键的是，高维矩阵的特征值排序更难，轨迹交叉更加频繁，贪心策略的失败率会显著上升。
常见做法是先对回比矩阵做一次平衡截断。平衡截断降维后的模型在低频段能保持原始系统的输入输出增益特性，特征轨迹的低频段形态基本保持一致，而高频段的差异可以通过频率加权修正。另一个技巧是使用复模态分解，把高频段的柔性问题解耦掉再分别处理。针对实际工程，如果只是课程设计验证，直接在原始模型上跑也够用；如果是高维工业系统，建议至少做一次最小实现降维，能省掉不少无用功。
5.3 频响求解加速：利用稀疏结构与并行池
对大规模稀疏状态空间模型，频响矩阵的求解是主要瓶颈。前面compute_loop_freq函数中用了左除，但如果A矩阵是稀疏的，并且你已知它的结构模式，可以先用分解一次性完成预处理，再在每个频点重复利用分解结果。MATLAB中可用分解对象封装LU分解，避免每个频点重复进行矩阵分解。
如果一台机器上有多个核心，还能用parfor代替for循环，把不同频点的频响计算和特征值分解分发给多个工作进程。但要注意parfor要求循环内的数据不能有依赖关系，在特征值排序部分不能直接套parfor，只能对频响计算部分做并行，排序和绘图仍然要串行完成。我实测过一个6阶模型，在8核并行下频响计算部分提速约5倍，但总体提速因为串行部分的存在只到了3倍左右，瓶颈反而转移到了绘图和排序阶段。
5.4 绘图刷新策略：避免800个频点逐点绘制的卡顿
直接把800个频点的散点全部绘制在一张图上，MATLAB的图形渲染会非常卡顿，尤其在高DPI屏幕和大量曲线叠加时。更合理的策略是先保存所有特征值轨迹的坐标数据，全部计算完成后用一次plot命令批量绘制，而不是在循环里逐点plot。
脚本里我采用的是每次plot一小段的策略，这是因为需要在循环中动态匹配颜色和连接线段。如果追求渲染速度，可以在循环内只做数据累积，退出循环后统一绘图。这个改动对代码整洁度有小幅影响，但对图形性能的提升非常明显——800个频点的2×2系统可能差别不大，但5×5系统800个频点就是4000个散点，逐点绘制的卡顿会让人崩溃。
6. 验证与进阶：闭环奇异值边界与频响残差检验
6.1 用闭环奇异值边界验证判据结论
广义奈奎斯特判据只回答稳定与否，没回答稳定裕度有多大。工程上更关心裕度，而系统的鲁棒稳定裕度和回比矩阵的奇异值边界直接相关。绘制广义奈奎斯特曲线后，我习惯叠加一组灵敏度函数S(s) = (I + L(s))^(-1)的最大奇异值曲线作为验证，如果最大奇异值在所有频段都小于某个阈值（比如2到3），说明闭环系统的输入输出增益被抑制在合理范围，特征轨迹判据得出的稳定性结论也是可信的。
实现方法很直接：沿用前面已经算好的Lfreq三维数组，逐频点计算inv(I + Lk)，再求svd取最大奇异值，最后用semilogx绘制以频率为横轴的曲线。如果在某个频段最大奇异值尖峰过高，说明该频段的反馈作用很弱，即使闭环稳定，系统的抗扰能力也差。这个验证并不复杂，代码大约十来行，资源里附带了完整的实现脚本。
6.2 检查频响残差以识别建模误差
资源里还有一个检测脚本用来评估模型不确定性对特征轨迹的影响。做法是在标称模型基础上叠加一个乘性不确定性Δ，重新绘制广义奈奎斯特曲线，然后对比标称轨迹和摄动轨迹的最大偏差。如果偏差过大，说明系统对模型误差敏感，这时候即使标称闭环稳定，实际系统也可能不稳定。
这个脚本适合在课程设计答辩或项目验收前做一次验证。把你自己系统的标称模型和实际工作点的辨识模型分别绘制曲线，叠加在同一张图上，观察特征轨迹穿越临界点的情况是否一致。在模拟项目X里我做过类似的验证，实测发现标称模型的高频段特征轨迹和摄动模型有明显偏离，后来就在控制器设计里加入了环路成形，把高频段增益压低，最终两条轨迹的一致性改善了很多。
6.3 特征轨迹法与其他稳定性判据的对比验证
想确认绘制结果没问题，最直接的办法是用其他方法交叉验证。对于低阶MIMO系统，可以直接计算det(I + L(s))的零点，通过求根函数求出闭环特征根，检查实部是否全为负，然后把结果与广义奈奎斯特判据的结论对比。两者一致说明判据应用正确，不一致就要怀疑特征值排序是否出错了。
另一种交叉验证是直接对闭环状态空间模型求极点。用ss(A - BKC, ...)构造闭环系统，然后计算闭环极点，如果所有实部都小于零，稳定性结论得到确认。这个方法虽然缺乏频域层面的直观性，但作为判断稳定性的底线是十分可靠的。从那以后我每次画完广义奈奎斯特曲线，都强制走一遍闭环极点验证，确认结论一致后才准备输出到文档或下一阶段的分析，省掉了很多后续无谓的调试，希望帮到你。
本文还有配套的精品资源，点击获取
