★ FEATURED ARTICLE 简介面向无线通信、信号处理方向的学习者与工程人员的MATLAB仿真素材聚焦多输入多输出系统的广义奈奎斯特曲线绘制用于观察开环传递函数矩阵在复平面上的稳定边界辅助判断MIMO系统的稳定性与频率、相位响应特征。压缩包共7个文件均为m脚本整体约5KB覆盖广义奈奎斯特绘制、符号传递函数转换、模型降阶与示例调用等环节便于逐段阅读源码并复现实验。已有3166人学习说明该主题在课程设计与科研入门中需求稳定。读者可据此理解nyqmimo函数的调用方式与参数含义结合示例修改系统矩阵、替换信道条件观察曲线变化并借助降阶与符号运算脚本排查建模问题形成从理论到MATLAB实现的完整练习路径适合具备控制或通信基础、希望动手验证MIMO稳定性分析的读者。1. 一个 2×2 回路单独调都稳一起投运就振荡某开发者调过一个 2×2 的过程对象两个回路分别用 PI 整定单回路的阶跃响应都很干净两个回路同时闭合之后输出开始以固定频率来回摆动而且频率不随设定值变化。翻单回路的 Bode 图幅值裕度 12 dB、相位裕度 45°怎么看都不该振。问题不在参数在判据本身经典 Nyquist 看的是单个标量回路对 -1 点的绕法MIMO 系统的稳定性由返回比矩阵 L(s)G(s)K(s) 的特征值轨迹决定任意一个通道的传递函数都不足以代表它。这份 MATLAB 资源干的就是这件事——把 L(jω) 在每个频点上的特征值算出来按频率连成特征轨迹画到复平面上再按广义奈奎斯特判据统计绕 -1 点的次数反推闭环右半平面极点个数。手上已经有 MIMO 模型、正在做分散控制或者多变量频域设计的开发者拿到就能对着自己的对象改模型跑。2. 项目骨架从模型入口到绘图管线的分层MATLAB 里画单回路 Nyquist 是一条nyquist(G)就完事MIMO 没有开箱即用的等价物因为「画什么」这件事本身要先定义清楚。这套代码的组织方式是把整条链条切成四段每段一个文件谁也不掺和谁的事。2.1 目录结构与职责划分mimo-nyquist/ ├── main_nyquist.m # 入口脚本模型 → 扫频 → 特征轨迹 → 判据 → 出图 ├── nyq/ # MATLAB 包目录调用时带 nyq. 前缀 │ ├── freqSweep.m # 频点生成、虚轴极点区域的采样处理 │ ├── retRatio.m # 返回比矩阵 L(jω) G(jω)K(jω) │ ├── charLoci.m # 逐频点特征分解 跨频点轨迹匹配 │ ├── winding.m # 绕数统计det(IL) 辐角法 │ └── plotLoci.m # 绘图、按频率着色、标注 -1 点 ├── models/ │ ├── model_tf.m # 传递函数矩阵入口 │ ├── model_ss.m # 状态空间入口 │ └── model_frd.m # 实测频响数据入口 ├── utils/ │ ├── scaleIO.m # 输入输出量纲归一化 │ └── closestPerm.m # 轨迹排序穷举或匈牙利匹配 └── tests/ └── test_cases.m # 已知结论的算例回归nyq是 MATLAB 的 package 机制目录名前面带加号里面的函数在外部要用nyq.charLoci(...)这样的全名调用。好处是命名空间隔离不会跟 Control System Toolbox 里已有的nyquist、pole、freqresp撞车——这一点在工具箱版本升级之后特别值钱。models/下面三个入口对应三种上游数据来源。这是整个项目里最该先想明白的一层抽象不管你的对象是解析建模的传递函数矩阵、状态空间实现还是从现场扫频仪导出的频响数据下游管线只认一种东西——形状是n×n×nw的复数三维数组第二个维度是频点数。2.2 三种模型入口统一成三维复数数组function L buildReturnRatio(G, K) % 把各种模型形式统一成返回比矩阵对象 % G: n×n 的 tf / ss / frdK: 可以是模型也可以是数值增益矩阵 if isnumeric(K) K tf(K); % 常数增益矩阵转成静态增益保证后面乘法合法 end L G * K; % 返回比矩阵MIMO 下乘法顺序不能反G*K ≠ K*G end function Ljw sampleReturnRatio(L, w) % 在频点 w 上取 L(jω)输出统一压成 n×n×nw 的三维数组 n size(L, 1); resp freqresp(L, w); if ndims(resp) 2 resp reshape(resp, n, n, 1); % 单个频点时 freqresp 会退化成二维 end Ljw reshape(resp, n, n, []); end第一段里isnumeric(K)那个分支别嫌啰嗦。实际调试时最常干的事就是改控制器增益矩阵比如把diag([1.2 0.8])换成diag([2.0 2.0])看轨迹怎么挪总不至于每次都tf(diag(...))。第二段里ndims(resp) 2的判断是血泪经验freqresp在只有一个频点时会返回n×n的二维矩阵而不是三维数组后面Ljw(:,:,k)直接索引报错报错信息还指向别处很容易看错方向。2.3 绘制管线的四步走整条管线是四步每一步的输入输出都是明确的改的时候可以只动其中一段。第一步生成频点。频率范围是个物理问题不是数值问题下限要覆盖最慢极点一个数量级以上上限要保证 L(jω) 已经衰减到接近零。第二步逐频点求 L(jω) 并做特征分解。第三步跨频点匹配特征值把属于同一条轨迹的点串起来。第四步统计绕数、定稳定性、出图。% main_nyquist.m 主干节选 w nyq.freqSweep(1e-2, 1e3, 1200); % 角频率 rad/s L buildReturnRatio(G, K); Ljw sampleReturnRatio(L, w); [lam, lamSorted] nyq.charLoci(Ljw); % 逐点求特征值 跨点匹配 [Ncw, phi] nyq.winding(Ljw, w); % 顺时针绕数 P sum(real(pole(L)) 0); % 开环右半平面极点数 Z P Ncw; % 闭环右半平面极点数 nyq.plotLoci(lamSorted, w, struct(Ncw, Ncw)); fprintf(P%d Ncw%d Z%d\n, P, Ncw, Z);freqresp走的是 Control System Toolbox 的模型求值对tf、ss、frd都能吃这也是把入口统一到这一层的原因。有实测频响数据的时候直接frd(mag, phase, w)构造后面的代码一行都不用改。管线里唯一跟“模型”打交道的就是sampleReturnRatio其余都只处理复数数组。3. 核心算法拆解从频响矩阵到特征轨迹这一章是整份资源最值得逐行读的部分。MIMO 的广义奈奎斯特判据听起来玄落地就是「每个频点算一次特征值串起来数圈」。3.1 逐频点特征分解eig 返回的顺序不可信eig对同一个矩阵求特征值返回顺序没有保证。L(jω₁)给的是[λa; λb]频率挪一格之后L(jω₂)可能给[λb; λa]。你要是直接按返回顺序把点连起来画出来的轨迹上会出现大跨度的横跳看上去像穿越了 -1 点实际上只是两条轨迹交换了位置。这是整个项目里最常见的翻车点。function [lam, lamSorted] charLoci(Ljw) % Ljw: n×n×nw 的复数三维数组 [~, n, nw] size(Ljw); lam zeros(n, nw); for k 1:nw lam(:, k) eig(Ljw(:, :, k)); % 每个频点单独求返回顺序是任意的 end lamSorted closestPerm(lam); % 相邻频点按距离做一一匹配 end循环里的eig单个频点调用就够了不必用批量接口因为真正的开销在后面的匹配上而匹配必须串行做每一列依赖前一列的结果。特征值的物理含义要强调一下λᵢ(jω)是返回比矩阵在各个方向上的“通道增益”复数有模有相。MIMO 系统的稳定性就由这些复数随频率在复平面上的走向决定。n 个特征值就有 n 条轨迹每条都要看。3.2 轨迹排序用相邻频点的最小总距离做匹配排序问题的本质是一个指派问题。n 阶一般不大2 到 6 阶最典型直接穷举 n! 种排列完全扛得住。function out closestPerm(lam) [n, nw] size(lam); out zeros(n, nw); out(:, 1) lam(:, 1); % 起点不定序后续以前一列为基准 P perms(1:n); % 全部排列n≤6 时最多 720 行 for k 2:nw prev out(:, k-1); cur lam(:, k); bestCost inf; bestPerm 1:n; for r 1:size(P, 1) idx P(r, :); cost sum(abs(cur(idx) - prev)); % 一一对应后的总复距离 if cost bestCost bestCost cost; bestPerm idx; end end out(:, k) cur(bestPerm); end endcur(idx)按排列重排当前频点的特征值跟上一列逐元素相减取模求和取总距离最小的那个排列。参数上唯一要留意的是n穷举是 n!n6 是 720n8 是 40320乘上频点数之后就会明显变慢。我自己实践里超过 6 阶就改用 MATLAB 自带的matchpairs把代价矩阵abs(cur - prev.)喂进去得到的结果是一样的复杂度低一个量级。matchpairs要求代价矩阵非负复数取模之后天然满足。还有一个容易忽略的细节这个匹配是贪心的、逐段独立的。某一段匹配错了后面会一路错下去。所以出图之后如果看到某条轨迹在某个频段突然折叠回来先回头查这一段附近的频率点是不是太稀。频率间隔太大两条轨迹在一格之内走过头交换了位置匹配算法也救不回来。3.3 绕数统计数 det(IL) 的相位不要一条条数圈这是这套代码里我最认可的一个设计选择。教科书上讲广义奈奎斯特判据说的是“n 条特征轨迹对 -1 点的净顺时针绕数之和”。直接照字面实现等于自己给自己挖坑轨迹从 -1 点附近擦过去的时候到底绕没绕数值上根本说不清angle在 ±π 处一抖就是一圈。更稳的等价做法是用辐角原理。特征值的性质保证det(I L) Πᵢ (1 λᵢ(jω))所以det(IL)的相位变化量等于所有(1λᵢ)相位变化量之和也就是所有特征轨迹绕 -1 点绕数的总和。把标量函数的绕数问题拿到标量上去算每次只处理一个复数序列unwrap就够了。function [Ncw, phi] winding(Ljw, w) % 用辐角原理统计 det(IL) 的绕数比逐条轨迹数圈稳得多 n size(Ljw, 1); nw numel(w); phi zeros(1, nw); for k 1:nw phi(k) angle(det(eye(n) Ljw(:, :, k))); end phi unwrap(phi); % 消除 ±2π 锯齿得到连续相位 dPhi phi(end) - phi(1); Ncw -round(dPhi / (2*pi)); % 逆时针为正取负得到顺时针绕数 end符号约定必须写死在注释里不然过两个月自己回来都看不懂。unwrap之后相位是连续增长的逆时针绕一圈相位涨 2π所以dPhi/(2π)是逆时针绕数广义奈奎斯特判据里 N 用的是顺时针取负号。判据本身是Z P N_cwZ 是闭环右半平面极点个数P 是开环右半平面极点个数闭环稳定等价于 Z 0。P sum(real(pole(L)) 0); Z P Ncw; stable (Z 0);pole(L)求的是返回比矩阵所有元素的公共极点对ss模型就是 A 矩阵的特征值。这一步有个前提L(s) 里不能有零极点相消否则pole会给出多余的开环极点P 算多了判据直接偏掉。还有一条适用范围要卡死L(s) 必须是方阵。输入输出个数不等的时候det(IL)都没有定义特征值也无从谈起。非方阵的 MIMO 稳定性得换奇异值或者结构奇异值那套工具这份代码不管。4. 绘图细节围线处理、着色与多回路对比算法算对了只是第一步图画错了一样会得出相反结论。这一章处理画图环节里几个真会咬人的地方。4.1 虚轴上的极点与频率采样Nyquist 围线是整条虚轴加上一个无穷大半圆。理论上如果 L(s) 在虚轴上有极点纯积分环节、未阻尼振荡环节都是围线必须往右做一个小半圆绕过去叫 indent。数值实现上不绕这一下freqresp在那个频率点直接返回Inf或者NaNdet出来是NaNangle出来还是NaNunwrap一碰到就整条废掉。function w sweepWithIndent(wlo, whi, np, wPole, r) % 生成扫频点wPole 是虚轴极点角频率列表r 是绕行半径 w logspace(log10(wlo), log10(whi), np); if nargin 4 || isempty(wPole), return; end mask true(size(w)); for p wPole(:) mask mask abs(w - p) r; % 把极点邻域整段挖掉 end w w(mask); fprintf(避开 %d 个虚轴极点剩余 %d 个频点\n, numel(wPole), numel(w)); end参数r的取值要跟极点位置和邻域频率点密度一起看一般是wPole的 1% 到 5%。太小了没避开太大了会把有效信息挖掉。说明白一点上面这段只是「挖点」严格意义上的 indent 还要在半圆上补采样把绕过极点的那段路径画出来。补采样的做法是把s沿j*p r*exp(1j*θ)θ 从 -π/2 走到 π/2逐个求G(s)*K(s)插到频率序列里对应位置。这套代码里freqSweep留了indent选项走的就是这条路默认关闭。如果你的对象有纯积分环节又拿不准把indent打开代价是频率序列变成非单调的绕数统计那段得改成按路径顺序累加相位而不是phi(end)-phi(1)。4.2 频率着色让「走向」在静态图里看得出来MIMO 特征轨迹跟单回路 Nyquist 最大的不同是单回路你只需要知道绕没绕特征轨迹你还想知道在哪一段频率绕的。所以按频率着色比用单一颜色画一条线信息量大得多。function plotLoci(lamSorted, w, opts) n size(lamSorted, 1); cmap parula(numel(w)); % 冷色低频暖色高频 hold on; for i 1:n for k 1:numel(w)-1 plot(real(lamSorted(i, k:k1)), imag(lamSorted(i, k:k1)), ... Color, cmap(k, :), LineWidth, 1.1); end plot(real(lamSorted(i,1)), imag(lamSorted(i,1)), o, ... MarkerEdgeColor, k, MarkerFaceColor, w); % 低频起点 end plot(-1, 0, rp, MarkerSize, 14, MarkerFaceColor, r); axis equal; grid on; box on; xlabel(\mathrm{Re}); ylabel(\mathrm{Im}); title(sprintf(特征轨迹 N_{cw}%d P%d Z%d, opts.Ncw, opts.P, opts.Z)); end逐段plot的做法可读性好代价是段数多的时候慢——1200 个频点乘 4 条轨迹就是 4800 次绘图调用出图能等上好几秒。真要快换成surface一次画完再用shading interp但代码就没这么好读了。日常调参时我会把频点降到 400 先用着定稿再跑 1200。axis equal不是可选项。复平面上实轴和虚轴比例不一致的话圆会画成椭圆绕没绕 -1 点用眼睛判断会出错。4.3 多组控制器对比与量纲问题分散控制最常见的需求是「换几组增益看看哪组把轨迹推得离 -1 点最远」。同一张图上叠几组曲线用线型区分比开几个 figure 来回切强。Ks {diag([1.2 0.8]), diag([2.0 2.0]), diag([0.6 1.5])}; styles {-, --, :}; figure(Color, w); hold on; for j 1:numel(Ks) L G * Ks{j}; Ljw sampleReturnRatio(L, w); lams closestPerm(realEig(Ljw)); [Ncw, ~] nyq.winding(Ljw, w); for i 1:size(lams, 1) plot(real(lams(i,:)), imag(lams(i,:)), styles{j}, LineWidth, 1.1); end fprintf(K%d: Ncw%d, Z%d\n, j, Ncw, sum(real(pole(L))0) Ncw); end量纲这件事得单独提一句。MIMO 里各通道的输入输出物理量常常不同一个是阀门开度 0-100%一个是转速 0-3000 rpm。直接把不同量纲的通道拼成一个矩阵求特征值得到的数在数学上没毛病在工程上没意义——特征值的大小被量纲最大的那个通道主导了。utils/scaleIO.m做的就是这件事给每个输入输出乘一个缩放因子把有效范围拉到同一量级再进算法。缩放因子不影响稳定性结论相似变换不改变特征值绕数但会影响你对图表可读性的判断。5. 避坑清单五个真会让人怀疑人生的现象5.1 现象与根因对照现象根本原因处理方式轨迹上出现大跨度横跳看着像穿越 -1eig返回顺序在频点间不一致相邻频点做最小距离匹配绕数算出来差 1改频点数又变了频率分辨率不足-1 附近相位跳变加密频点用det(IL)辐角法低频段一片空白判据给“稳定”扫频下限高于最慢极点下限取最慢极点 1/10NaN/Inf混进结果整张图空白虚轴极点处freqresp发散极点邻域挖点或 indent 绕行手算结论和代码结论差一个符号顺时针/逆时针约定没统一注释写死约定用闭环极点校验5.2 逐条排查与自检脚本第一条轨迹跳变。现象是图上某两条曲线在一个频点处突然交叉换位走过去再换回来。原因是eig顺序问题不是物理现象。解决办法就是 3.2 节的匹配但匹配完还得回头看一眼最大单步跳变量超过阈值就说明频率点太稀了匹配算法也拿不准。第二条绕数抖动。现象是同一个系统频点从 800 加到 1200Ncw从 -1 变成 0。原因几乎都是轨迹从 -1 点附近很近的地方掠过采样点刚好落在两侧。用det(IL)的相位连续变化去算能缓解但不能根治。真正的解法是在|1det(IL)|最小的那段频率区间里插值加密甚至可以把最小距离打印出来看看。[~, idx] min(abs(arrayfun((k) det(eye(n)Ljw(:,:,k)), 1:numel(w)))); fprintf(最接近 -1 的频点: w%.4g, 距离%.4g\n, ... w(idx), abs(det(eye(n)Ljw(:,:,idx))));第三条频率范围不够。这是新手最容易犯的。默认logspace(-2, 3)听着挺宽但对象里要是有个时间常数 1000 秒的慢模态1e-2 rad/s 完全没覆盖到。反过来高频端如果没衰减到零det(IL)的相位还在缓慢漂移绕数就会算进一个不真实的增量。所以高频收敛要断言。tailErr abs(det(eye(n) Ljw(:, :, end)) - 1); assert(tailErr 0.05, ... 高频段 det(IL) 未收敛到 1当前 %.3g扩大扫频上限, tailErr);第四条NaN和Inf。带积分环节的对象一定有这个问题。第一反应是w(w ~ 0)把零点剔掉这没用因为极点不一定在 0。得用pole(L)找出所有虚轴上abs(real(p)) 1e-8的极点频率列表里把它们的邻域整段跳过或者按 4.1 节走 indent。跳过之后记得在图上标一下断点不然看的人以为轨迹就那么断的。第五条符号约定。这份代码里定的是Ncw -round(dPhi/(2*pi))顺时针为正判据Z P Ncw。换个来源的教材可能用的是逆时针为正、Z P - N看起来只是符号差别结果会稳定判反。防这个的办法只有一个每次改完绕数那段代码拿一个已知稳定的算例比如L diag(1/(s1), 1/(s2))P0Z 应该也是 0跑一遍输出对不上就是符号搞反了。别靠脑子记。6. 二次开发把脚本做成能进回归测试的判稳函数原项目入口是脚本clear all开头那种改模型得进去改几行。自己用没问题一旦要拿它去批量检查一组工况或者接进自动化流程就必须函数化。改造的核心是三件事输入参数化、输出结构化、结论可交叉验证。第一件事用arguments块把接口卡死。模型、控制器、频率范围、收敛容差都从参数进函数体里不出现任何硬编码的数字。function res mimoNyquist(G, K, opt) % MIMO 广义奈奎斯特判稳 % 输出 res: Ncw(顺时针绕数) / P / Z / stable / lam / w / tailErr arguments G K eye(size(G, 1)) opt.Wlo (1,1) double 1e-2 opt.Whi (1,1) double 1e3 opt.Np (1,1) double 1200 opt.Tol (1,1) double 0.05 end n size(G, 1); w nyq.freqSweep(opt.Wlo, opt.Whi, opt.Np, ... imag(pole(G*K)) 1e-8, 0.02*opt.Whi); Ljw sampleReturnRatio(buildReturnRatio(G, K), w); [~, lamSorted] nyq.charLoci(Ljw); [Ncw, ~] nyq.winding(Ljw, w); P sum(real(pole(G*K)) 0); Z P Ncw; tailErr abs(det(eye(n) Ljw(:,:,end)) - 1); res struct(Ncw, Ncw, P, P, Z, Z, stable, Z 0, ... lam, lamSorted, w, w, tailErr, tailErr); endarguments块是 R2019b 之后才有的老版本换成inputParser或者手写nargin分支逻辑一样。第二件事也是最关键的一件事拿闭环极点当后悔药。绕数法算出来的 Z本质上是在间接回答「闭环有没有右半平面极点」那就直接把闭环极点求出来对一下对不上说明前面某一步有问题。function checkAgainstPoles(G, K) % 用闭环极点交叉验证绕数结论跑回归测试时必走 res mimoNyquist(G, K); L G * K; clp pole(feedback(L, eye(size(L,1)))); % 默认负反馈 Zref sum(real(clp) 1e-8); assert(Zref res.Z, ... 绕数法给出 Z%d闭环极点实际 Z%d检查频率范围与绕数符号, ... res.Z, Zref); fprintf(交叉验证通过Z%d高频收敛误差 %.2e\n, res.Z, res.tailErr); endfeedback(L, eye(n))默认负反馈如果你的回路是正反馈接法要加第三个参数1这个参数位置很多人记反。real(clp) 1e-8用容差而不是严格大于零是因为数值误差会让本应在虚轴上的极点飘出一点点。第三件事回归用例。我一般准备三类开环稳定且闭环稳定P0, Z0、开环不稳定但闭环能镇住P1, Z0绕数必须是 -1、故意调坏的增益Z0。三组都过了才敢拿它去跑真实工况。维度原始脚本函数化之后模型切换改代码第 8 行传参批量扫描控制器增益手动重跑arrayfun一行结论可信度靠眼睛看图assert交叉验证接自动化流程做不到返回结构体可序列化最后一个常用技巧arrayfun扫增益。把控制器增益做成网格一次跑完把使Z0的那些点筛出来画成参数平面上的稳定域——这是分散控制里最实用的图比一张张看特征轨迹效率高得多。g1 linspace(0.2, 3.0, 40); g2 linspace(0.2, 3.0, 40); [G1, G2] meshgrid(g1, g2); Zmap arrayfun((a, b) mimoNyquist(G, diag([a b])).Z, G1, G2); contourf(G1, G2, Zmap, LineColor, none); colorbar; xlabel(k_1); ylabel(k_2); title(闭环右半平面极点个数);有一次我图省事直接拿脚本改了两行就交出去图上看着稳稳的结果现场一投运就振回头查是logspace下限写成了 1 而不是 1e-2慢模态整个被漏在扫频窗口外面绕数自然算成 0。从那以后我每次改完频率范围都强制先跑一遍闭环极点交叉验证再出图对不上就不看曲线——图会骗人极点不会。希望这套流程能帮到你少走一次这种弯路。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
