MIMO广义奈奎斯特曲线绘制与特征轨迹判稳实践

简介面向无线通信、信号处理方向的学习者与研究人员这份资源聚焦多输入多输出MIMO系统的广义奈奎斯特曲线绘制用MATLAB脚本演示如何由开环传递函数矩阵计算稳定边界并在复平面上呈现多输入多输出之间的频率响应与相位关系帮助读者判断系统稳定性、理解信道容量分析中的关键判据。压缩包内共7个文件均为m脚本涵盖MIMO广义奈奎斯特绘制、符号传递函数与状态空间转换、模型降阶及示例调用等模块整体约5KB体积轻量但环节完整便于逐步阅读源码并结合理论对照曲线含义。资源已有3166人学习说明其在课程实验与课题验证中有一定参考价值。对于需要快速搭建MIMO稳定性分析流程、复用绘图脚本或补充模型转换工具的读者可直接在MATLAB中运行示例修改系统参数或模拟不同信道环境观察曲线变化从而加深对广义奈奎斯特判据与MIMO系统性能特点的理解。1. 从一次“单回路都稳、整机却自激”说起MIMO 广义奈奎斯特曲线到底画的是什么调试多输入多输出MIMO回路时最容易出错的一幕是把每个回路单独做阶跃超调、上升时间、相位裕度都漂亮一合上交叉通道输出就开始等幅振荡。原因不复杂——单回路裕度只看对角元耦合项被当成了不存在。要判断这类系统的闭环稳定性得把整个回路矩阵 L(s)G(s)K(s) 当作一个整体看它的特征值 λi(jω) 在复平面上扫出的轨迹离 -1 点多远、净绕了几圈这就是广义奈奎斯特曲线也叫特征轨迹。用 MATLAB 把这条曲线画出来本身不难难的是画对、把绕数数对、让曲线不断线这几处基本都得自己写。2. 判据与技术选型特征轨迹为什么替代不了逐回路裕度2.1 广义奈奎斯特判据到底在判什么单回路那套 Nyquist 判据的推导起点是开环传递函数 1L(s)多变量系统里对应的量是 det(IL(s))。把行列式展开成特征值形式det(IL) ∏(1λi)其中 λi 是 L(jω) 的特征值。于是闭环特征方程的零点、开环 RHP 极点、以及 det(IL(jω)) 对原点的绕数之间仍然满足辐角原理那套关系只不过“对 -1 点的绕数”要按每条特征轨迹分别统计再求和闭环稳定 ⟺ Σ Ni P其中 Ni 是第 i 条特征轨迹对 -1 点的净绕数逆时针为正P 是 L(s) 在右半平面的极点数重极点按重数计。这条判据有两个容易被忽略的前提。一是 L(s) 必须是方阵即输入数等于输出数非方阵系统得先做 squared-down 或者直接对 det(IL) 求绕数但后者只能给出稳定与否给不出“离不稳定还有多远”的几何直觉。二是 P 要按 L(s) 的极点算如果 G 和 K 之间存在零极点对消而 det(IL) 里这部分被约掉了P 的取值就会算错判据直接失效——我见过一次模型是从辨识结果里直接拼的分子分母有几乎相同的因子P 算出 2实际闭环是稳的图形上也看不出异常排查了大半天。真正让工程上愿意用特征轨迹的原因是它把“稳定裕度”变成了看得见的几何量每条曲线离 -1 点的最近距离 dmin、曲线穿过单位圆时的频率、曲线在 -1 点附近的走向都能直接读出来。当某条轨迹贴着 -1 点掠过时你就知道这个方向上的增益再抬一点就会出事而不像逐回路裕度那样四个回路都显示 6 dB 裕度、合起来却只剩 1 dB。2.2 MATLAB 里三条可选路线与取舍在 MATLAB 里做多变量频域分析能走的路不止一条差别在于你要的是图、是数、还是裕度指标。路线关键函数输出适用场景逐点特征分解freqrespeig特征轨迹坐标、可自定义配对与着色出图、数绕数、教学演示行列式奈奎斯特手算det(eye(n)L(jw))后画复平面轨迹只有总绕数快速验证判据、写单测盘裕度 / 结构化奇异值diskmargin、mu、musyn同时允许增益与相位摄动的裕度鲁棒性评估、和轨迹图交叉验证第一条路线是这次的重点因为只有它能拿到每条轨迹的原始坐标后面做配对、着色、自动判据全都依赖这些中间量。第二条路线的价值在于交叉验证det(IL) 的绕数应当等于各特征轨迹绕数之和两个数对不上基本可以断定配对逻辑或者频率网格出了问题。第三条路线属于“结论层”diskmargin会返回一个盘型裕度它考虑的是所有回路同时、同方向地摄动增益和相位这个指标比逐回路裕度保守也更贴近实际但它是黑匣子式的输出出问题时没法告诉你问题出在哪个频率、哪个方向。这里必须点破一个高频误解直接对 MIMO 系统调用nyquist(L)MATLAB 画出来的是每个输入输出对的元素级 Nyquist 曲线ny×nu 张图不是特征轨迹。这两者名字像、长得也像但判据完全不同——元素级曲线对 -1 的绕数在多变量系统里没有稳定性含义。我带的实习生第一次做这个就是拿元素级图去解释自激现象结论全反了。3. 工程骨架一套可复用的绘图工具的目录与数据流3.1 目录结构一次性的脚本能出图但改两次参数之后就没人敢动了。我习惯把它拆成一个带包目录的小工程函数名统一挂在mnyq命名空间下调用时写mnyq.charLoci(...)避免和控制工具箱里的函数撞名。mimo_nyq/ ├── main_demo.m # 一键跑通建模 - 算轨迹 - 画图 - 打印绕数 ├── mnyq/ # 包目录所有算法函数挂在这个命名空间下 │ ├── charLoci.m # L(jw) 逐点特征分解返回特征值与特征向量 │ ├── pairLoci.m # 特征值配对保证曲线跨频点连续 │ ├── windingNumber.m # 绕数统计 到 -1 点的最小距离 │ ├── plotLoci.m # 频率着色绘图、临界点与单位圆标注 │ └── gainScan.m # 增益扫描输出稳定域表格 ├── models/ │ ├── model2x2.m # 2x2 强耦合对象用来验证算法正确性 │ └── model4x4.m # 4x4 过程对象用来看高维下的糊图问题 ├── config/ │ └── plotCfg.m # 频率范围、点数、配色、坐标范围等集中配置 └── tests/ └── test_pairing.m # 人造交叉特征值专门验证配对逻辑把模型单独抽成函数是因为判据验证需要反复换对象把配置集中到一个 cfg 结构体是因为频率范围和点数这两个参数几乎每次调试都要改散落在脚本里迟早会改漏一处。tests/test_pairing.m看起来有点小题大做但配对逻辑是全工程最容易出隐性错误的地方——它错了图还是能画出来只是曲线偷偷换了道肉眼很难发现。3.2 从模型到曲线的六步数据流整条链路的中间变量形状变化很容易记混尤其是freqresp返回的三维数组索引顺序踩过一次坑就会记一辈子。步骤变量维度说明建模G、Kn×n 的tf/ss建议最终转成状态空间数值更稳组成回路Ln×nL G*K注意是矩阵乘不是逐元素频率网格wnw×1logspace生成单位 rad/s频率响应frn×n×nwfreqresp(L,w)的返回顺序特征分解lam、Vnw×n、n×n×nw每个频点一次eig配对后lamnw×n列索引代表“第几条轨迹”判据与绘图N、dmin1×n每条轨迹的绕数与最近距离数据流里唯一有分支的是最后一步要图就走plotLoci要结论就走windingNumber两个函数吃的是同一份配对后的lam。这种设计的好处是图上的几何直觉和判据的数值结论永远来自同一批数据不会出现“图看着稳、代码判不稳”的对不上账。4. 核心代码拆解从 L(jω) 到能看懂的特征轨迹4.1 建模与频率网格先把网格铺对% 2x2 被控对象 对角控制器故意保留较强的交叉通道 s tf(s); G [ 1/(s1), 0.6/(s2); 0.4/(s3), 2/(s0.5) ]; K diag([1.5, 0.8]); % 先只看比例控制耦合影响最明显 L G * K; % 回路传递矩阵注意是矩阵乘法 L minreal(L, 1e-6); % 消去数值上的零极点对避免响应毛刺 % 频率网格对数均匀铺底 在 0.1~10 rad/s 区间加密 w [ logspace(-3, -1, 80), linspace(0.1, 10, 400), logspace(1, 4, 120) ]; w unique(w(:)); % 去重并转成列向量minreal那一步不是可选项。辨识或者手工拼接出来的模型分子分母常带 1e-8 量级的近似公因子freqresp在少数频点上会算出明显偏离的值反映到特征轨迹上就是曲线上突然冒出一个尖刺看图的人会以为系统在某个频率有谐振。容差 1e-6 是经验值太小了消不掉太大了会把真实的相近极点也吃掉。频率网格用了三段拼接而不是单个logspace。原因是特征轨迹穿越 -1 点附近的区域通常集中在穿越频率周围对数均匀采样在这段区间点数偏少绕数统计容易漏掉一次小回环。低频段用 80 点是因为低频特征值变化平缓不需要密高频段用logspace(1,4,120)是因为曲线在高频会收敛到常值矩阵的特征值点太密只是浪费计算。真要省事logspace(-3,4,2000)也能用代价是eig被调用两千次4x4 以上的系统会明显卡顿。4.2 逐点特征分解freqresp 返回的三维数组怎么用function [lam, V, w] charLoci(L, w) %CHARLOCI 计算 L(jw) 的特征值特征轨迹与对应特征向量 % lam : nw x n 复数矩阵第 k 行是频率 w(k) 处的 n 个特征值 % V : n x n x nw第 k 页的列向量是 lam(k,:) 对应的特征向量 n size(L, 1); nw numel(w); fr freqresp(L, w); % 返回 n x n x nw 的三维数组 lam zeros(nw, n); V zeros(n, n, nw); for k 1:nw [Vk, D] eig(fr(:,:,k)); % 单频点特征分解 lam(k,:) diag(D).; % 转置成行方便按频率堆叠 V(:,:,k) Vk; end endfreqresp的返回顺序是“输出×输入×频率”第三个维度才是频率这一点和bode返回的mag维度一致但和很多人的直觉相反。写成fr(:,:,k)取第 k 个频点的 n×n 复矩阵直接丢给eig即可。必须同时保留特征向量V这是后面配对的前提。eig返回的特征值顺序是算法内部决定的同一个系统在相邻两个频点上第 1 个特征值可能属于完全不同的分支——eig本身不保证任何连续性。如果只留特征值画出来的曲线会在某些频点突然“换道”看起来像是断成了几截实际上数据没错是索引错了。顺带说一个数值细节eig对非对称实矩阵返回的是复数特征值理论上 L(jω) 在实频率上的特征值不保证共轭成对出现除非 L 是实系数矩阵在实频率取值此时才共轭成对。所以判断“某条轨迹是否穿到右半平面”不能简单看虚部符号得看实部。4.3 特征值配对让曲线不断线的关键一步function [lam, V] pairLoci(lam, V) %PAIRLOCI 按特征向量相关性把相邻频点的特征值接上消除顺序跳变 % 思路v_i(jw_k) 与 v_j(jw_{k-1}) 的内积模越大越可能是同一条分支 [nw, n] size(lam); for k 2:nw C abs(V(:,:,k) * V(:,:,k-1)); % n x nC(i,j)|v_i^k, v_j^{k-1}| usedR false(n,1); usedC false(n,1); perm 1:n; for t 1:n Cc C; Cc(usedR,:) -inf; Cc(:,usedC) -inf; % 已占用的行列屏蔽 [~, idx] max(Cc(:)); [r, c] ind2sub(size(Cc), idx); % r: 新频点序号, c: 旧频点槽位 perm(c) r; % 槽位 c 放第 r 个特征值 usedR(r) true; usedC(c) true; end lam(k,:) lam(k, perm); V(:,:,k) V(:, perm, k); end end这是整套代码里最值得单独写一个函数的地方。核心假设是频率变化足够小时同一条特征轨迹的特征向量方向也只会小幅旋转所以用相邻频点特征向量的内积模作为“相似度”做一次贪心指派就能把顺序对齐。用内积模而不是实部是因为特征向量带任意相位v和e^{jθ}v表示同一个方向取模消掉这个自由度。Cc(usedR,:) -inf这几行是在做互斥一旦某个新序号被分配出去它所在的整行和对应槽位的整列都要屏蔽否则贪心会重复指派。贪心不保证全局最优遇到两条轨迹在某个频点几乎重合特征值简并时可能选错但代价是错一次只影响局部一段不会像排序法那样整条曲线重排。要更稳可以换成匈牙利算法matchpairs(-C, 1e3)R2019a 之后才有内部实现是全局最优指派代价是每步 O(n³)n 大于 8 之后比贪心慢一个数量级。配对做完之后还有个隐性收益V的顺序也一起对齐了后面要算特征向量的相位变化率、或者做模态分析直接就能用。4.4 绕数与临界点把“看图”变成“判据”function [N, dmin] windingNumber(lam, w, target) %WINDINGNUMBER 统计每条特征轨迹绕 target默认 -1的净圈数 % N : 1 x n逆时针为正、顺时针为负的净绕数 % dmin : 所有轨迹到 target 的最小距离衡量稳定裕度 if nargin 3, target -1; end [nw, n] size(lam); N zeros(1, n); dmin inf; for i 1:n z lam(:,i) - target; % 平移到以临界点为原点 d abs(z); dmin min(dmin, min(d)); dth diff(unwrap(angle(z))); % 解卷绕后的相位增量 N(i) round(sum(dth) / (2*pi)); % 净绕数 总相位增量 / 2pi end endunwrap是这里的命门。直接对angle的结果做差分角度在 ±π 处会跳变 2π累加出来的总相位增量就完全不对了。unwrap把跳变补回去sum(dth)/(2π)才是真实的净绕数round用来吸收浮点累积误差正常情况下结果应该非常接近整数如果算出来是 0.5 这种值说明频率网格太粗轨迹在某个频段绕了半圈就被跳过去了。还有一个绕数统计的范围问题上面这段只处理了 ω 从 0 到 ∞ 的半条轨迹。完整的 Nyquist 围线包含负频率部分而 λ(-jω) 与 λ(jω) 共轭所以完整绕数等于半条轨迹绕数的两倍再加上无穷远处的圆弧贡献。对严格真有理的 L(s)分母阶次高于分子ω→∞ 时 L→0圆弧收缩成一个点不贡献绕数对双真系统L(∞) 是个常值矩阵那一段得按常值矩阵的特征值单独补一段路径。代码里没写这段是因为示例对象都是严格真的换对象时必须自己确认这一点不然判据会差一圈。dmin这个输出用起来比绕数更顺手。绕数只有 0、±1、±2 几个离散值而dmin是连续的比值小于 0.3 就说明曲线贴着临界点走哪怕当前还是稳定的增益稍微一动就会翻车。4.5 绘图频率着色与坐标等比例function plotLoci(lam, w, cfg) %PLOTLOCI 以频率着色的方式画特征轨迹并标注临界点与单位圆 figure(Color,w); hold on; cmap parula(numel(w)); % 新版可换 turbo对比度更高 for i 1:size(lam,2) for k 1:numel(w)-1 % 逐段着色才能体现频率方向 plot(real(lam(k:k1,i)), imag(lam(k:k1,i)), -, ... Color, cmap(k,:), LineWidth, 1.4); end plot(real(lam(1,i)), imag(lam(1,i)), o, MarkerSize,5); % 起点 plot(real(lam(end,i)), imag(lam(end,i)), s, MarkerSize,5); % 高频端 end plot(-1, 0, k, MarkerSize, 14, LineWidth, 2); % 临界点 th linspace(0, 2*pi, 400); plot(cos(th), sin(th), k:); % 单位圆参考 axis equal; grid on; xlabel(Re); ylabel(Im); caxis(log10([w(1) w(end)])); colorbar; title(sprintf(广义奈奎斯特曲线%d 条特征轨迹, size(lam,2))); endaxis equal是必须的。默认坐标轴会按数据范围自动拉伸一条本来离 -1 点还有 0.5 距离的曲线在横向被压扁之后看起来像是直接穿过去了看图的人会得出完全相反的结论。加上axis equal实轴虚轴同比例几何距离才有意义。逐段plot而不是一次性画整条是为了让颜色随频率变化——对数频率下曲线的疏密本身就携带信息低频段点稀、高频段点密颜色一叠上去哪一段是穿越区一眼就能认出来。缺点是段数多的时候渲染慢2000 个频点、4 条轨迹就是 8000 次plot调用出图能等上好几秒要批量出图就换成surface或者patch一次性铺牺牲一点可控性换速度。caxis在较新的版本里改名叫clim老代码搬过去会报错改成clim(log10([w(1) w(end)]))即可。5. 避坑与排查五条踩过的坑5.1 现象曲线莫名其妙断成几截或者两条轨迹在某个频率交换了位置原因是eig不保证相邻频点的特征值顺序一致而排序法按实部或按相角排序在特征值接近简并时会判错两条轨迹整段换道图上看起来就像断了。解决办法就是 4.3 节的向量配对用特征向量相似度做指派。判断有没有配对成功有个简单办法把lam每列相邻点的差值取模如果某处出现远大于平均值的跳变比如比中位数大两个数量级说明那里没接上。tests/test_pairing.m就是干这个的构造一对交叉的特征值函数跑一遍看配对结果和解析解是否一致。5.2 现象绕数算出来不是整数或者和 det(IL) 的绕数对不上两个原因。一是频率网格不够密轨迹在穿越区绕的小回环被采样的直线段“切掉”了二是漏算了负频率和无穷远弧段的贡献。处理办法分两步先用det(eye(n)L(jw))单独算一遍总绕数作为真值再拿特征轨迹的绕数之和去比对不上就加密网格特别是穿越频率附近用linspace局部补点。我在 4x4 对象上遇到过差一圈的情况最后发现是高频段logspace(1,4,120)太稀一条轨迹在 500 rad/s 附近有个很小的回环被完全跳过补到 400 点之后才对上。5.3 现象图上四条轨迹挤成一团完全看不出哪条离 -1 最近高维系统常见。n 大于 3 之后特征轨迹在复平面上的分布会非常拥挤尤其是对象各通道增益差异大的时候有的轨迹贴近原点有的跑到几十以外坐标范围一拉开临界点附近全糊了。解决办法是先做通道缩放把 L 的输入输出按稳态增益归一化让各通道量级接近再画图实在挤不开就分面画subplot每张图只放一两条轨迹配上各自的dmin标注。缩放的意义不只是好看——广义奈奎斯特判据本身对相似变换不变缩放不改变绕数所以可以放心做。5.4 现象把nyquist(L)的图当成广义奈奎斯特曲线这是最典型的一个。nyquist(L)对 MIMO 系统画的是元素级曲线2x2 系统出四张子图每张对应一个输入输出对判据和特征轨迹完全不是一回事。区分方法很直接元素级图的曲线数量是 ny×nu特征轨迹的数量恒等于 n再不行就看代码里有没有eig(freqresp(...))这一步没有就不是。写报告的时候我一般会两张图都放元素级用来说明耦合强度特征轨迹用来下稳定性结论读者不会混淆。5.5 现象模型改了一个参数判据从稳定变成不稳定但阶跃响应看着正常多半是零极点对消导致 P 算错。L G*K之后如果没做minreal分子分母里残留的近似公因子会让pole(L)返回一堆虚假的右半平面极点P 从 0 变成 2判据立刻翻转。排查顺序先zpk(L)看一眼零极点分布有没有相距小于 1e-4 的对子再做minreal(L, 1e-6)最后用闭环极点max(real(pole(feedback(L, eye(n)))))交叉验证。闭环极点全在左半平面、而判据说不行那一定是 P 算错了不要反过来怀疑判据。提示判据代码和闭环极点验证最好永远成对出现。绕数只是个数字闭环极点才是最终事实两者对不上时优先信极点然后回头查频率网格、配对和 P 的计算。6. 二次开发从一张图到一个能批量判稳的小工具单张图解决的是“看一眼”的需求。实际项目里更常见的场景是控制器结构定了要扫一遍增益看稳定边界在哪或者对象有多个工况每个工况都要判一次。这时候把 4.4 节的绕数函数往外扩一层就行。function T gainScan(G, K0, gains, w) %GAINSCAN 扫描标量增益 k用广义奈奎斯特判据判断闭环稳定性 n size(G,1); T table(); for k gains L minreal(G * (k*K0), 1e-6); [lam, V] mnyq.charLoci(L, w); [lam, ~] mnyq.pairLoci(lam, V); N mnyq.windingNumber(lam, w, -1); P sum(real(pole(L)) 0); % 开环 RHP 极点数 Ntot sum(N); % 各轨迹绕数之和 Tcl feedback(L, eye(n)); poleOK all(real(pole(Tcl)) -1e-8); % 闭环极点真值 T [T; table(k, P, Ntot, Ntot P, poleOK)]; end end关键在最后两列Ntot P是判据给出的结论poleOK是闭环极点的真实结论两列并排放在表里跑十次只要有一次不一致立刻就能定位到是哪个增益点出了问题。我一般还会在表里加一列min(abs(lam(:)1))记录该增益下离临界点的最近距离稳定域边界附近这一列会趋近于零比单纯看 true/false 有信息量得多。进一步可以做的扩展有三个方向。一是把w做成自适应网格先粗扫一遍找出dmin的极小值所在频段再在那个频段局部加密重算能省掉一大半计算量。二是和diskmargin交叉验证两者的结论趋势应当一致diskmargin偏保守如果出现盘裕度说稳定而绕数说不稳定那基本可以确定是 P 或网格的问题。三是出报告把特征轨迹图、绕数表、稳定边界图拼成一张用exportgraphics直接导出 PDF评审的时候比一堆散图省事。有个细节值得单独说gainScan里每次循环都重建了L对tf对象来说这步很慢。如果G和K0都是状态空间模型把G和K0提前转成ss、用k直接缩放K0的 B、C 矩阵循环速度能快一个数量级4x4 对象扫 50 个增益点tf版本要跑两分多钟状态空间版本十几秒就完事。踩过最疼的一次是扫增益时忘了在循环里做minreal某个增益点恰好让两个极点靠得极近pole判定翻车判据给出的稳定域边界偏了整整 0.4。从那以后我每次写这类批量扫描都强制在循环体内做一次模型化简并且在结果表里保留闭环极点验证那一列——宁可多算一遍也不想再拿着一张边界位置错误的图去开会。希望这些能帮到你少走一遍我走过的弯路。本文还有配套的精品资源点击获取