简介面向无线通信、信号处理方向的学习者与研究人员这份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(IL(s)) 直接相关。广义奈奎斯特判据的表述是令 s 沿奈奎斯特围线虚轴加右半平面大圆弧走一圈L(s) 的特征值 λ_i(s) 在复平面上画出若干条轨迹这些轨迹对临界点 (-1, 0) 的净包围圈数 N满足 Z N P其中 P 是 L(s) 在右半平面的极点数Z 是闭环右半平面极点数。Z0 即闭环稳定。这里有两个容易翻车的细节。第一包围圈数要按“净”算逆时针为正、顺时针为负方向搞反结论直接颠倒。第二P 必须从开环模型里数出来不能默认是 0——很多实际系统开环就带右半平面极点比如某些欠阻尼对象或不稳定平台P 不为零时判据的结论会完全不同。资源里如果只画图不给 P 和 N 的计数那图只能看个热闹。单回路判据是这套东西在 L(s) 退化为标量时的特例此时特征值就是 L(s) 本身轨迹只有一条。MIMO 下特征值有多条每条都要数包围且特征值轨迹在特征值重合点附近会出现“交换”现象排序错乱会让曲线看起来断裂——这是绘图层面最常见的视觉坑后面避坑章节会专门讲。2.2 为什么用特征值而不是行列式理论上 det(IL(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/10wmax 取最高关注频率的 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*(s1)); G12 0.8/(s2); G21 0.5/(s3); G22 6/((s0.5)*(s2)); G [G11 G12; G21 G22]; % 控制器示例对角比例积分 K diag([1.2*(11/(0.5*s)), 0.9*(11/(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/s3000 点 [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 点是经验值若曲线在穿越区仍有折角加到 5000w的下限不要低于系统最低极点实部的 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即 -10junwrap对扫频点密度敏感点太稀时相位跳变超过 π 会误判这也是为什么前面强调 3000 点起步。最后Z N PZ 为 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 处理导致结论错误现象开环明明有右半平面极点判据却按 P0 算得出“稳定”的错误结论。原因直接套 ZN 而没数开环极点或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×nevalfr返回矩阵后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)*(11/(0.5*s)), gains(b)*(11/(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、NZmap0标出稳定区。参数上扫参步长 0.1 是精度和耗时的折中2000 点扫频在 25×25 网格下大约几十秒量级若嫌慢可先粗扫再在边界附近细扫。这张稳定域图比单张奈奎斯特曲线有用得多——它直接告诉你两个通道增益的可行组合边界而单回路裕度给不出这个信息。验证方法上我习惯用闭环极点做交叉校验对稳定域内的点用feedback(L, eye(2))求闭环极点确认全部实部为负对边界点闭环极点应恰好贴近虚轴。这一步是后悔药——广义奈奎斯特的计数一旦有符号或密度问题闭环极点会立刻暴露比盯着曲线数圈数可靠。另外特征值轨迹在重特征值附近对参数极敏感扫参时若发现某网格点 Z 突变先怀疑是数值问题而不是系统真不稳定把该点扫频加密重算一遍再下结论。从那以后我每次交付频域稳定性结论都强制走一遍“特征轨迹图 P/N 显式计数 闭环极点交叉校验”三件套少一步都不敢签字。希望帮到你。本文还有配套的精品资源点击获取