简介面向无线通信、信号处理方向的学习者与研究人员这份MATLAB代码资源聚焦多输入多输出MIMO系统的广义奈奎斯特曲线绘制帮助评估系统稳定性与频率响应特性。压缩包共7个文件全部为.m脚本整体约5KB涵盖传递函数符号化转换、模型降阶、广义奈奎斯特绘图及示例主程序等模块便于按流程理解从建模到成图的完整链路。已有3163人学习下载说明其在教学与工程分析中具备一定参考价值。读者可通过运行示例脚本观察复数传递函数矩阵在复平面上的稳定边界表现并结合参数调整与信道环境模拟加深对MIMO系统相位特性与稳定条件的理解同时提升MATLAB编程与通信理论结合应用的能力。1. 从一张“绕圈”的曲线说起MIMO 广义奈奎斯特到底在画什么做多变量反馈设计的同行大概都遇到过这种场面单回路里看奈奎斯特曲线一眼就能判断稳定裕度可一旦系统变成两输入两输出把每个通道的奈奎斯特图分别画出来曲线看着都离临界点很远闭环却偏偏发散。问题不在图画错了而在于多变量系统的稳定性不能靠“逐个通道看”来判断——回路之间存在耦合一个通道的增益变化会牵动另一个通道的等效开环特性。广义奈奎斯特曲线就是为解决这件事而生的它把 MIMO 开环传递矩阵 L(s) 在每个频率点的特征值 λ_i(jω) 取出来把这些特征值轨迹画在同一张复平面上用它们绕 (-1, 0) 点的圈数来判断闭环稳定性。这份 MATLAB 资源围绕的正是这条曲线的绘制流程适合正在做多变量频域设计、需要把稳定性判据落到图上的控制方向开发者。2. 技术底座与目录结构这套 MATLAB 工程由哪些文件撑起来2.1 为什么选频域特征值轨迹而不是时域仿真时域仿真能告诉你“这一组参数下发散”但说不清“离发散还有多远、往哪个方向调参数更安全”。广义奈奎斯特判据的价值在于给出裕度信息特征值轨迹离临界点的最近距离、绕圈方向、穿越负实轴的频率区间这些都是可以直接指导参数整定的量。MATLAB 在这个场景下的优势是eig、freqresp、nyquistplot这些函数现成传递矩阵可以直接用tf或ss对象拼出来不用自己写频响计算。常见做法是把被控对象写成状态空间形式控制器写成传递矩阵两者相乘得到开环 L(s)再在关心的频段上扫点求特征值。需要提前说清一个概念边界广义奈奎斯特判据针对的是方阵开环传递矩阵输入输出维数必须相等。如果你的系统是 2 入 3 出这类非方情形得先做输入或输出变换把它方阵化否则特征值无从谈起。这一点在动手写代码前就要确认不然后面所有图都是错的。2.2 工程目录的典型组织方式这类绘制工程一般不会只有一个脚本按职责拆开更利于二次开发。下面是我拆这类项目时常见的目录骨架文件名按功能命名你可以对照自己手上的资源包核对mimo_nyquist/ ├── main_nyquist.m % 主入口定义对象、控制器调用绘制流程 ├── model/ │ ├── build_plant.m % 构造被控对象传递矩阵 G(s) │ └── build_controller.m % 构造控制器矩阵 K(s) ├── core/ │ ├── open_loop.m % 计算开环 L(s)G(s)K(s) │ ├── eig_sweep.m % 频段扫描 逐点特征值分解 │ └── winding_number.m % 绕圈数统计与临界点距离 ├── plot/ │ ├── plot_gnq.m % 广义奈奎斯特曲线绘制 │ └── plot_margin.m % 裕度标注最近距离、穿越频率 └── config/ └── freq_range.m % 频率范围与采样密度配置这个划分的核心思路是“模型构造 / 数值计算 / 绘图”三层分离。模型层只负责把 G(s)、K(s) 拼出来计算层只吃传递矩阵吐特征值序列绘图层只吃数据画图。好处是换一个被控对象时只改model/下的两个文件绘图和绕圈数统计完全复用反过来想换绘图风格也不碰数值代码。我一般会把频率范围和采样点数单独抽到config/因为这两个参数对曲线形状影响极大集中管理比散落在脚本里好排查。2.3 关键函数依赖与版本注意整个流程依赖 Control System Toolboxfreqresp、tf、ss、eig都在这个工具箱里。如果你的 MATLAB 没装这个工具箱tf和freqresp会直接报未定义函数这是新手最容易卡住的第一步。另外nyquistplot在较新版本里对 MIMO 对象的支持有变化自己用plot画特征值轨迹反而更可控——这也是为什么这类工程通常不直接调nyquistplot而是手动算特征值再plot。3. 核心计算链路从开环矩阵到特征值轨迹的逐段实现3.1 构造开环传递矩阵 L(s)一切从开环开始。设被控对象 G(s) 是 2×2控制器 K(s) 也是 2×2开环 L(s) G(s)K(s)。用 MATLAB 的传递矩阵写法每个元素是一个tf对象拼成矩阵后矩阵乘法自动按传递函数运算% build_plant.m —— 构造 2x2 被控对象传递矩阵 s tf(s); % 对角通道两个一阶惯性环节时间常数不同 g11 1 / (s 1); g22 2 / (0.5*s 1); % 耦合通道交叉项体现回路间相互影响 g12 0.3 / (s 2); g21 0.2 / (s 3); G [g11, g12; g21, g22]; % 控制器简单的比例 积分对角结构 k11 1.5 0.8/s; k22 1.2 0.5/s; K [k11, 0; 0, k22]; % 开环传递矩阵 L G * K;这段代码的关键点有三个。第一s tf(s)声明了传递函数变量后面所有1/(s1)这类写法才会被识别为传递函数而不是数值除法漏了这一行会得到一堆标量而不是tf对象。第二交叉项g12、g21不能省它们正是“单通道看着稳、整体却发散”的根源去掉耦合项再画广义奈奎斯特你会得到一个过于乐观的结论。第三控制器用对角结构是为了先验证基础流程实际工程里 K 可能是满阵矩阵乘法G*K照样成立不需要改代码结构。参数上时间常数和增益决定了特征值轨迹的走向。g11的极点 -1 和g22的极点 -2 会让低频段特征值偏向实轴正方向耦合项系数 0.3、0.2 越大两条特征值轨迹分得越开。调这几个数你能直观看到曲线形状怎么变这是理解判据最快的办法。3.2 频段扫描与逐点特征值分解有了 L(s)下一步是在对数分布的频率点上求每个点的特征值。频率范围要覆盖系统的主要动态下限低于最慢极点一个数量级上限高于最快极点一个数量级否则曲线两端会“缺头少尾”绕圈数统计算不准。% eig_sweep.m —— 频段扫描并逐点求特征值 function [lambda_all, w] eig_sweep(L, wmin, wmax, npts) % 对数分布频率点避免低频段采样过疏 w logspace(log10(wmin), log10(wmax), npts); % 计算开环频率响应返回 2x2xnpts 的三维数组 Lresp freqresp(L, w); nout size(Lresp, 1); lambda_all zeros(nout, numel(w)); % 每列存一个频率点的全部特征值 for k 1:numel(w) Lk Lresp(:, :, k); % 取出第 k 个频率点的复矩阵 lam eig(Lk); % 特征值分解 % 按实部排序保证相邻频率点的特征值轨迹连续 [~, idx] sort(real(lam)); lambda_all(:, k) lam(idx); end end这里有个血泪经验eig返回的特征值顺序是任意的不同频率点之间顺序可能跳变直接画出来会得到一堆交叉乱线根本看不出轨迹。解决办法是按实部或按与上一频率点特征值的最近距离排序让同一支轨迹在相邻频率点上“接得上”。上面用实部排序是最简单的做法对大多数工程对象够用如果两支轨迹实部很接近、排序仍会跳就改成“每个特征值匹配上一列中距离最近的那个”代码多几行但更稳。freqresp返回的是三维数组第一维输出、第二维输入、第三维频率点取Lresp(:,:,k)才是第 k 个频率点的复矩阵。这个索引顺序搞反是常见翻车点症状是eig报维度不匹配或者算出来的特征值完全不对。npts一般取 500 到 2000点太少曲线折线感重、绕圈数统计会漏点太多绘图卡顿但精度提升有限1000 是个稳妥的起点。3.3 绕圈数统计与临界点距离广义奈奎斯特判据的核心是数特征值轨迹绕 (-1, 0) 的净圈数。数值上不能靠“看”得用辐角累积法把每个频率点相对临界点的复向量角度算出来相邻点角度差做解缠绕累加后除以 2π 就是绕圈数。% winding_number.m —— 统计单支特征值轨迹绕 (-1,0) 的净圈数 function N winding_number(lam_track, crit) if nargin 2 crit -1 0i; % 默认临界点 -1j0 end % 相对临界点的复向量 rel lam_track - crit; % 相邻点角度差unwrap 处理 ±pi 跳变 dtheta diff(unwrap(angle(rel))); % 净绕圈数 总角度变化 / 2π N sum(dtheta) / (2*pi); endunwrap这一步不能省。angle返回的是 (-π, π] 区间内的值轨迹穿过负实轴时角度会从 π 跳到 -π不 unwrap 就会把一次正常穿越算成反向绕圈圈数直接错。crit参数留成可配置是因为有些教材把临界点定义在 (-1, 0)有些在原点附近做变换留个口子方便对照。除了圈数实际调参更关心“离临界点多近”。在eig_sweep的结果上再算一个最小距离% 所有特征值轨迹到临界点的最小距离作为稳定裕度的粗略指标 all_lam lambda_all(:); min_dist min(abs(all_lam - (-1 0i)));这个距离不是严格的稳定裕度定义多变量系统有更严谨的 μ 分析但作为快速筛查足够用距离越小参数越接近危险区。我一般会在绘图时把这个最近点和对应频率标出来调参时盯着它走。4. 绘图与结果解读把特征值轨迹画成能读的图4.1 广义奈奎斯特曲线绘制数据齐了绘图本身不复杂难点在“让图能读”多支轨迹要区分颜色、临界点要标出、单位圆和负实轴要画参考线。% plot_gnq.m —— 绘制广义奈奎斯特曲线 function plot_gnq(lambda_all, w) figure; hold on; grid on; ntrk size(lambda_all, 1); colors lines(ntrk); % 每支轨迹一种颜色 for i 1:ntrk plot(real(lambda_all(i,:)), imag(lambda_all(i,:)), ... Color, colors(i,:), LineWidth, 1.2, ... DisplayName, sprintf(\\lambda_%d, i)); % 标出低频起点方便判断轨迹走向 plot(real(lambda_all(i,1)), imag(lambda_all(i,1)), ... o, Color, colors(i,:), MarkerFaceColor, colors(i,:)); end % 临界点 plot(-1, 0, rx, MarkerSize, 12, LineWidth, 2, ... DisplayName, critical point (-1,0)); % 单位圆参考 theta linspace(0, 2*pi, 200); plot(cos(theta), sin(theta), k--, LineWidth, 0.5); axis equal; xlabel(Real); ylabel(Imag); title(Generalized Nyquist Plot (eigenvalue loci)); legend(Location, best); endaxis equal必须加否则实轴虚轴比例不一致圆看起来是椭圆绕圈关系会被视觉扭曲。低频起点用实心圆点标出来是因为轨迹方向随频率增大往哪走对判断绕圈方向很关键光看一条线分不清头尾。lines(ntrk)生成区分度足够的颜色序列轨迹支数多时比手动指定颜色省事。4.2 结果怎么读三个判断动作图出来之后按顺序做三件事。第一看每支轨迹是否包围 (-1, 0)包围则对应闭环不稳定在开环无右半平面极点的前提下这是判据的直接结论。第二看轨迹离临界点的最近距离和对应频率距离小说明该频段裕度不足调参优先动这个频段的增益。第三看轨迹穿越负实轴的位置穿越点越靠左越接近 -1相位裕度越紧张。需要提醒的是广义奈奎斯特判据的完整形式还要考虑开环右半平面极点数 P闭环不稳定极点数 Z P NN 为绕圈数符号约定因教材而异。如果你的对象本身开环不稳定比如倒立摆类P 不为零光看“包不包围”会得出错误结论。这类工程里通常会在winding_number旁边再放一个pole(L)统计右半平面极点数两个数一起看才完整。4.3 和单回路奈奎斯特图的对照把同一系统的单通道奈奎斯特图和广义奈奎斯特图并排看是理解耦合影响最直观的方式。单通道图里每个对角元素的曲线可能都远离 -1但广义图里特征值轨迹却绕了进去——这个反差正是多变量系统“通道裕度不等于系统裕度”的证据。做二次开发时我建议保留一个对照绘图函数把两种图放一个 figure 的两个 subplot 里汇报和自查都用得上。5. 避坑与排查这类绘制工程最容易翻车的五个地方5.1 曲线画出来是乱麻轨迹互相穿插现象绘图结果是一团交叉的线看不出任何连续轨迹。原因eig返回顺序在相邻频率点间跳变没有做排序或匹配。解决在eig_sweep里按实部排序或者用“最近邻匹配上一列特征值”的方式重排排序后仍乱检查频率点是否过疏导致轨迹跳变过大把npts提到 2000 再试。5.2 绕圈数算出来是小数或明显不对现象winding_number返回 0.5、1.7 这类非整数。原因angle的 ±π 跳变没做unwrap或者频率范围没覆盖完整轨迹导致首尾不闭合。解决确认unwrap在diff之前调用把频率下限再降一个数量级、上限再升一个数量级让轨迹两端都收敛到稳定值检查crit是否和判据约定一致。5.3 非方阵系统直接报错现象eig报“矩阵必须为方阵”。原因开环 L(s) 是 m×n 且 m≠n特征值无定义。解决广义奈奎斯特判据只适用于方阵开环非方系统要先做输入/输出变换比如取部分输出做反馈、或引入权重矩阵做方阵化这一步没有通用公式得按具体系统结构设计别硬套。5.4 换了对象后图完全不对但代码没报错现象曲线形状离谱和预期动态对不上。原因freqresp返回的三维数组索引顺序搞反把输入维当成了频率维或者tf变量s在函数内部重新声明导致作用域混乱。解决在eig_sweep里打印size(Lresp)确认是[nout, nin, nfreq]把s tf(s)统一放在模型构造函数内部不要跨函数共享。5.5 低频段曲线发散到无穷远现象轨迹在低频端冲向无穷。原因控制器含积分环节1/s低频增益趋于无穷特征值模值随之爆炸。解决这是物理事实不是 bug绘图时把频率下限设在积分起作用之上一点或者对幅值做截断显示判据分析时低频段本来就要单独处理别指望一张图看全。6. 二次开发把绘制流程改造成自己的稳定性检查工具拿到这套代码后最值得做的改造不是换绘图配色而是把它变成一个能批量跑参数的检查函数。我的做法是把main_nyquist.m里的流程包成一个函数输入是对象参数和控制器参数输出是绕圈数、最小距离、是否稳定三个量这样就能套for循环扫参数网格找出稳定域边界。% check_stability.m —— 批量稳定性检查的封装 function [is_stable, min_dist, N_total] check_stability(p) % p 是参数结构体增益、时间常数、耦合系数 G build_plant(p); K build_controller(p); L G * K; [lam, ~] eig_sweep(L, 1e-2, 1e3, 1000); % 统计所有轨迹的净绕圈数之和 N_total 0; for i 1:size(lam, 1) N_total N_total winding_number(lam(i,:)); end % 开环右半平面极点数 P sum(real(pole(L)) 0); % 闭环稳定判据Z P N_totalZ 为 0 则稳定 is_stable (P round(N_total)) 0; min_dist min(abs(lam(:) - (-1 0i))); end这个封装的关键是把“看图判断”变成“数值判断”round(N_total)处理浮点误差P从pole(L)直接统计。有了它扫参数就是两层循环的事把is_stable为真的参数点收集起来稳定域边界自然浮现。验证方法也简单拿一组已知稳定的参数跑一遍is_stable应为真、min_dist应明显大于零再故意把增益调大到发散is_stable应翻假。两个方向都对封装才算可信。还有一个容易忽略的改造点频率范围不该写死。不同对象的动态时间尺度差几个数量级1e-2到1e3只适合秒级对象。我一般从pole(G)和zero(G)的模值自动推频率范围取最小极点模值的百分之一到最大极点模值的百倍这样换对象不用手调。这个自动推范围的逻辑放进config/freq_range.m和手动配置并存默认走自动、特殊情况手动覆盖。从那以后我每次把这类频域判据代码交给别人用之前都强制拿一个已知答案的简单系统跑一遍对照——单输入单输出退化成经典奈奎斯特结果必须和手算一致这个回归测试能挡掉大部分索引和符号约定的低级错误。希望帮到你。本文还有配套的精品资源点击获取