简介：面向多输入多输出通信系统稳定性分析与教学验证的MATLAB代码包，适合无线通信、信号处理方向的研究生、教师及工程技术人员使用。内容围绕广义奈奎斯特曲线绘制展开，涉及开环传递函数矩阵建模、符号推导与模型降阶等环节，帮助读者在复平面上判断系统稳定边界与相位特性。压缩包共7个文件，均为.m脚本，体积约5KB，涵盖主绘图函数、示例脚本及配套的符号转换与降阶辅助函数，结构紧凑，便于逐行阅读与二次修改。目前已有3166人学习，说明其在相关教学与科研场景中具备一定参考价值。读者可直接运行示例观察不同参数与信道条件下的曲线变化，并对照源码理解各步骤的数学含义，从而把理论公式、MATLAB编程与稳定性判据串联起来，形成可复用的分析流程。

1. 三个回路各自都稳，合起来却发散了：MIMO 稳定性判据缺的那块拼图

三个各自稳定的单回路，接进同一台多输入多输出设备之后整体振荡——这种场景在 matlab 里能复现，在台架上也能复现。直觉上每条回路的开环奈奎斯特曲线都离 −1 点远得很，闭环能差到哪去？问题恰恰出在「单独」两个字上：MIMO 系统的稳定性由整条特征轨迹共同决定，断开其余回路去测其中一条，测到的只是「其他回路还开着」的时候那个局部情况。多输入多输出本质是 m×m 的传递矩阵，耦合项让每个输出同时受多个输入影响，相位裕度和增益裕度就不再是标量，而变成一组随频率旋转的数。

想真正判断闭环稳不稳，得把 m 条开环特征轨迹一并放到复平面上，观察它们是否作为整体绕过了 −1 点。这套判据叫广义奈奎斯特判据，围绕它画出来的曲线就是广义奈奎斯特曲线。这份 matlab 资源要解决的，正是「怎么把这条曲线画对、并且从图上读出一个可信的整数结论」这件事。

2. 广义奈奎斯特判据：为什么盯住 det(I+L) 而不是逐个回路

2.1 闭环极点藏在 det(I+L(s)) 的零点里

单位负反馈下，开环传递矩阵记作 $L(s)=G(s)K(s)$，$G$ 是被控对象，$K$ 是控制器，两者都是 $m\times m$ 的有理函数矩阵。闭环传递矩阵是 $(I+L)^{-1}L$，它的极点只可能来自两处：$(I+L)^{-1}$ 的奇点，以及 $L$ 自身零极点对消后残留的部分。前者由 $\det(I+L(s))=0$ 决定，后者在无对消假设下不贡献右半平面极点。

于是闭环稳定性就被压缩成一个标量有理函数的零点位置问题：$\det(I+L(s))$ 在右半平面有没有零点。这一步转换是整个判据的地基，也是它和单回路奈奎斯特的分水岭。单回路里判的是 $1+L(s)$ 这条曲线绕不绕原点，多回路里判的是 $\det(I+L(s))$ 这条曲线绕不绕原点，形式上一模一样，含义从「一条回路」变成了「整个矩阵的联合行为」。

把 $\det(I+L(s))$ 写成 $N_c(s)/D_L(s)$，$D_L$ 是 $L$ 各元素分母的最小公倍式，$N_c$ 就是闭环特征多项式（相差一个常数因子）。开环不稳定极点个数 $P$ 等于 $D_L$ 在右半平面的根数，闭环不稳定极点个数 $Z$ 等于 $N_c$ 在右半平面的根数。对 $\det(I+L)$ 沿右半平面围线用幅角原理，得到 $Z-P=N$，$N$ 是曲线顺时针绕原点的圈数。

约定符号必须一开始就固定，否则后面全乱。工程上习惯取「顺时针为正」，因为顺时针绕圈对应右半平面零点多于极点。这个约定一旦写进代码注释，就不要在中途改。

2.2 乘积结构：绕 −1 与绕原点是同一件事

$\det(I+L)$ 有一个非常好用的展开：$\det(I+L)=\prod_{i=1}^{m}(1+\lambda_i(L))$，$\lambda_i$ 是 $L$ 的特征值。每个因子 $1+\lambda_i$ 在复平面上就是从 $-1$ 指向 $\lambda_i$ 的向量。所以 $\det(I+L)$ 的相位等于各因子相位之和，模等于各因子模之积。

这个乘积结构直接给出两种等价的画法。第一种，直接画 $\det(I+L(j\omega))$ 这一条曲线，数它绕原点的圈数，判据是 $Z=P+N_{\det}$。第二种，把 $m$ 条 $\lambda_i(L(j\omega))$ 轨迹画出来，数它们各自绕 $-1$ 点的圈数之和，判据是 $Z=P+\sum_i N_i$。两条路在数学上完全等价，因为「某个 $1+\lambda_i$ 转过原点」和「某个 $\lambda_i$ 转过 $-1$」是同一件事。

等价不等于等效。行列式法的优点是只要一条曲线，视觉上干净；缺点是行列式是连乘，任意一个因子接近零都会把整条曲线拽到原点附近，数值相对误差被放大，而且曲线绕圈时你无法判断是哪个模态干的。特征值法的优点是能定位到具体模态，哪个回路、哪个方向在逼近 $-1$ 一目了然，调控制器时这个信息很值钱；缺点是特征值本身在重根附近不可微，轨迹会出现真正的交叉，连线逻辑要额外处理。

我一般两个都画：特征值轨迹用来看「谁在惹事」，行列式曲线用来做最终判定，两者绕数对不上就说明采样或数值出了问题。

2.3 P 从哪来：被控对象自身的右半平面极点

$P$ 是 $\det(I+L)$ 在右半平面的极点数，在无零极点对消的前提下等于 $L$ 的右半平面极点数。$L=GK$，所以 $P$ 由被控对象和控制器共同决定。绝大多数人被这个量误导过：以为单位反馈、控制器稳定，$P$ 就是零。

只有被控对象本身稳定、控制器也稳定时 $P=0$ 才成立。一旦对象存在右半平面极点——某些机械结构、某些化工过程的线性化模型都有——$P$ 就非零，此时哪怕曲线绕了 −1 一圈，系统照样可能稳定，因为 $Z=P+N$ 里 $P$ 把它抵消掉了。反过来，$P=0$ 而曲线顺时针绕了一圈，那 $Z=1$，闭环必有一个右半平面极点。

实际计算 $P$ 的稳妥做法是直接求 $L(s)$ 所有元素分母的根，统计实部大于零的个数，注意去掉重根。用状态空间模型更省事：如果 $G$ 有状态矩阵 $A$，$P$ 就是 $A$ 的特征值中实部为正的个数，前提是 $G$ 的实现是最小实现。非最小实现里那些不可控不可观的模态会被算进 $P$，导致判据整体偏移。

提示：$P$ 算错的后果比采样不够密严重得多，因为它是一个整数的整体平移。先确认 $P$，再谈曲线。

3. 扫频与特征值配对：奈奎斯特曲线能不能连对

3.1 对数网格为什么会在谐振峰处骗你

画频率响应，习惯上取 $\omega$ 的对数等间距网格，从 $10^{-2}$ 到 $10^{3}$ 取几百个点。这对 Bode 图够用，对奈奎斯特曲线经常不够。原因是曲线在 $-1$ 附近的曲率由 $\sigma_{\min}(I+L)$ 的极小值决定，而在谐振峰附近，特征值随频率的变化率可以大到相邻两个频点之间的相位增量超过 $\pi$。

一旦单步相位增量超过 $\pi$，相位累加法就会把方向判反，绕数直接错一个整数。更隐蔽的情况是增量接近 $\pi$，方向没错但曲线在两步之间实际绕过了 $-1$，图上看着没绕，判定结果却是绕了。

实用的加密策略是两段式。先按对数网格粗扫一遍，算出每个频点上的 $d(\omega)=\min_i|1+\lambda_i(L(j\omega))|$，找 $d$ 的局部极小值，在极小值附近用几何中点二分插入新频点，直到相邻点的相位增量都小于某个阈值（比如 $\pi/4$）。这比全频段均匀加密省一个数量级的计算量，因为绝大多数频段曲线走得很平稳。

function w = refine_grid(Lfun, w0, thresh) % Lfun : 输入频率向量，返回 [m m nw] 的回差矩阵采样 % w0 : 初始对数网格 % thresh: 允许的最大单步相位增量（弧度） w = w0; for iter = 1:20 % 最多加密 20 轮，防止死循环 L = Lfun(w); lam = eigen_trajectories(L); % 已做频率方向配对 dphi = abs(diff(unwrap(angle(lam.'), [], 1), 1, 1)); bad = max(dphi, [], 2) > thresh; % 每个频段的最坏相位增量 if ~any(bad), break; end wm = sqrt(w(bad) .* w([false; bad]));% 几何中点，保持对数分布 w = sort([w, wm]); end end

这段代码的关键在第 8 行：unwrap沿频率方向展开相位，diff求相邻步长，取每条轨迹的最大值。bad标记的是「这一段不可信」，而不是「这一段错了」。几何中点用sqrt(w1*w2)是因为对数网格上中点在几何意义下才有意义，算术中点会把点堆到高频端。

3.2 eig 返回的顺序没有语义，配对得自己做

这是画广义奈奎斯特曲线最容易翻车的地方。eig内部走的是 Schur 分解，返回特征值的顺序由算法决定，同一族矩阵在相邻频率上完全可能把第 1 个和第 3 个特征值互换。直接plot(real(lam(1,:)), imag(lam(1,:)))画出来的曲线会在某几个频点上出现横跨复平面的竖线，看起来像曲线自己绕了一大圈——绕数立刻错。

正确的做法是在频率方向做轨迹配对，把每个频点的特征值分配给上一频点中最接近的那一条轨迹。

function lam = eigen_trajectories(L) % L : [m m nw] 回差矩阵采样，nw 为频点数 % lam : [m nw]，每一行是一条沿频率连续的轨迹 [m, ~, nw] = size(L); lam = zeros(m, nw); lam(:,1) = sort(eig(L(:,:,1))); % 首点排序只为可复现 for k = 2:nw cur = eig(L(:,:,k)); prev = lam(:,k-1); C = abs(prev - cur.'); % 代价矩阵 C(i,j)=|prev_i - cur_j| lam(:,k) = cur(greedy_assign(C)); end end function M = greedy_assign(C) % 最小代价一一分配，返回 M(i) 表示第 i 条轨迹接哪个新特征值 [m, ~] = size(C); M = zeros(m,1); used = false(1,m); for i = 1:m Ci = C(i,:); Ci(used) = inf; % 已占用的列禁止再选 [~, j] = min(Ci); M(i) = j; used(j) = true; end end

C(i,j)是上一频点第 $i$ 条轨迹到当前频点第 $j$ 个特征值的欧氏距离，greedy_assign做的是按行贪心的一一匹配。严格意义上这应该用匈牙利算法求全局最优，但频率步长足够小时贪心解和最优解几乎总是一致，代价小很多。Ci(used) = inf那一行是防止两条轨迹抢同一个特征值，丢掉它会在重根附近产生重复分配。

这里有个容易绕进去的概念：即使配对错了，绕数在数学上也不该变。因为 $\lambda$ 集合的并集才是决定绕数的对象，分支如何连接不影响整体是否包围 $-1$——前提是轨迹不穿过 $-1$。所以配对的作用主要是两件事：让图能看，让基于单条轨迹的算法（比如逐条算相位）不出错。真正做稳定性判定时，用行列式法或对全部轨迹求和可以绕开配对问题。理解到这层，就知道配对是工程手段，不是数学必需。

3.3 用最小奇异值给曲线可信度打分

特征值在重根附近不可微，数值上表现为两列特征向量接近共线，eig的结果对扰动极其敏感。判断某个频点的特征值能不能信，看 $\sigma_{\min}(I+L(j\omega))$：它是 $I+L$ 到奇异矩阵的最小距离，也近似是特征值 $\lambda_i$ 到 $-1$ 的最小距离下界（在 $L$ 正规时才严格相等，一般矩阵上有差距，但趋势可用）。

function [smin, kappa] = cred_score(L) % L : [m m nw] 回差矩阵采样 % smin : 每个频点的最小奇异值，越小越接近奇异，越不可信 % kappa : 条件数，越大说明特征值方向越不稳定 nw = size(L,3); smin = zeros(1,nw); kappa = zeros(1,nw); for k = 1:nw s = svd(L(:,:,k) + eye(size(L,1))); smin(k) = s(end); kappa(k) = s(1) / max(s(end), eps); end end

smin接近零的频点要单独标记出来，这些点上det(I+L)的相对误差会被放大到不可接受，行列式曲线的位置不可信。经验阈值是smin < 1e-3 * smax时放弃行列式法，改用逐轨迹相位法并把采样加密一个量级。kappa用来判断矩阵本身是不是接近亏损，条件数超过 $10^{8}$ 时eig给出的特征向量基本没有意义，但特征值作为集合仍然可用——这一点常被误解成「条件数大就不能算特征值」，实际上不能算的是特征向量和单条轨迹的归属。

4. 绕数怎么数：把图上的「看着像」变成可判定的整数

4.1 相位累加法与 (−π, π] 折算

最直观的绕数算法是累加相位。曲线 $z(\omega)$ 相对参考点 $z_{\text{ref}}$ 的相位是 $\angle(z-z_{\text{ref}})$，沿频率逐点累加相邻相位差，折算到 $(-\pi,\pi]$ 再求和，除以 $2\pi$ 就是绕数。

function [N, dmax] = winding_phase(z, ref) % z : 复数折线，已按频率递增排好序 % ref : 参考点，行列式法取 0，特征值法取 -1 % N : 逆时针为正的绕数（正频率半段，未乘镜像系数） % dmax : 相邻段最大相位增量，用于判断采样可信度 zz = z(:) - ref; d = angle(zz(2:end)) - angle(zz(1:end-1)); d = mod(d + pi, 2*pi) - pi; % 折算到 (-π, π] N = sum(d) / (2*pi); dmax = max(abs(d)); end

第 6 行的mod(d+pi,2*pi)-pi是全部精度的关键。直接对angle的差值求和，在相位跨过 $\pm\pi$ 时会跳出 $2\pi$ 的跳变，导致绕数偏大。先加 $\pi$、对 $2\pi$ 取模、再减 $\pi$，等价于把差值折到最短路径上，跨象限时不会出错。

dmax是自检量。它接近 $\pi$ 说明某一段的采样不足，绕数不可信；它恰好等于 $\pi$ 说明曲线正好穿过参考点，此时绕数无定义，系统处在临界稳定边界上。这两个信息比绕数本身更早暴露问题。

注意符号约定：winding_phase返回逆时针为正，而 2.1 里的判据用的是顺时针为正。最终判定时取负号，或者把判据写成 $Z = P - N$。这种符号在代码里一定要写注释，我一个同事曾经因为符号搞反，把一个本来稳定的系统判成不稳定，白改了三天控制器。

4.2 有向交点法：不依赖相位展开的备选

相位法有一个隐患：unwrap类逻辑对噪声敏感。如果曲线本身有高频抖动（比如控制器里有微分环节），相位增量会在真值附近震荡，累加误差可能吃掉小半个整数。备选方案是有向交点法：从参考点向正实轴方向发一条射线，数折线穿过这条射线的次数，按穿越方向定正负。

function N = winding_cross(z, ref) % 有向交点法：射线从 ref 指向 +∞ 实轴方向 zz = z(:) - ref; x = real(zz); y = imag(zz); N = 0; for k = 1:numel(x)-1 y1 = y(k); y2 = y(k+1); if (y1 > 0) ~= (y2 > 0) % 两端点在实轴两侧 xi = x(k) - y1 * (x(k+1)-x(k)) / (y2-y1); % 线性插值求交点 if xi > 0 % 交点必须在正半轴 N = N + sign(y2 - y1); % 向上穿 +1，向下穿 -1 end end end N = -N; % 对齐「逆时针为正」约定 end

第 5 行的(y1 > 0) ~= (y2 > 0)判断两端点是否分居实轴两侧，这个写法把「端点在实轴上」的情况排除在外，避免重复计数。第 7 行的线性插值假设折线段足够短、可以近似为直线，这也是交点法对采样密度敏感的地方。第 8 行xi > 0是射线法而非整条实轴法的关键，只统计正方向那一侧。

交点法的好处是只依赖每一段的几何位置，不做跨段的相位累加，误差不累积。坏处是它需要曲线是闭合的——正频率段的曲线要和镜像段拼起来才能算总绕数。如果只算正频率半段，交点法给出的值和相位法差一个系数，不能混用。

4.3 自适应加密与临界穿越的判定

绕数算法跑出结果之后，必须再用dmax做一次可信度检查。一旦dmax超过阈值，就在超限频段插入新点重算，直到收敛。

function [N, ok] = winding_trusted(Lfun, w0, ref, thresh) % 循环加密直到相位增量全部低于阈值，返回可信绕数 w = w0; ok = false; N = 0; for iter = 1:25 L = Lfun(w); lam = eigen_trajectories(L); [N, dmax] = winding_phase(lam(:), ref); % 全部轨迹拼成一或多条分别算 if dmax < thresh, ok = true; break; end wm = sqrt(w(1:end-1) .* w(2:end)); w = sort([w, wm]); end end

这段代码里thresh一般取 $\pi/4$ 到 $\pi/6$。加密 25 轮还不收敛，通常意味着曲线正好穿过参考点，也就是系统处在临界状态，此时绕数为半整数或无定义，任何取整都是错的。这种情况要停下来看 $\sigma_{\min}(I+L)$ 是否接近零，而不是继续加压网格。

还有一个频率端点的问题。当 $L(s)$ 含积分环节时，$s=0$ 是极点，$\det(I+L(j0))$ 发散，正频率段不能从 0 开始扫，得从 $10^{-4}$ 之类的下限起。当 $\det(I+L(j\omega))$ 在 $\omega\to\infty$ 不趋于常数时，高频端要做同样的截断。这两个截断不影响统计 RHP 包围次数的整数值——根据小圆弧或大圆弧的渐近方向，包围次数要么不变，要么只差一个可以解析确定的项。但很多画出来的曲线在极限处直接跳到无穷再回来，肉眼看不出来这一跳对应多少次包围；遇到这种情况，宁可回到状态空间模型用精确的闭环极点核对，也不要硬数图形。

关于镜像，正频率段的绕数要乘 2 才是完整的 RHP 围线绕数。因为实系数有理函数满足 $L(-j\omega)=\overline{L(j\omega)}$，负频率段是正频率段的共轭反向，两者贡献相同，相加就是加倍。这一步非常容易被漏掉，漏掉之后 $Z=P+N$ 会把「不稳定」和「稳定」判反，而且反得很有规律——你可能会发现模型明明发散，判据却说稳定。

5. 排查：曲线画出来像那么回事，结论却是反的

5.1 图形层面：曲线连得对，眼睛才信得过

现象一：曲线在某一段出现横跨复平面的直线，看起来像突然绕了一大圈。

原因：eig在相邻频点返回的顺序发生交换，而代码里没做轨迹配对，直接按特征值下标连线。重根附近这个问题几乎必然出现。

解决：在频率方向做最近邻配对（见 3.2 的eigen_trajectories），或者在绘图前用sort按实部排序——后者只是把突跳从随机位置挪到固定位置，不能真正修好，只能用于快速看一眼。正式判定必须做配对。

现象二：图上明显没绕过 −1，判定结果却说绕了一圈。

原因：频率网格在谐振峰附近跨过了整段圆弧。曲线在两个采样点之间实际绕了 −1 一整圈又回来，采样点上完全看不出来。

解决：用|1+\lambda_i|的局部极小值定位关键频段，在极小值附近用几何中点加密，直到单步相位增量低于 $\pi/4$。更省事的办法是看smin曲线，smin的谷底就是必须加密的地方。

现象三：曲线在 $-1$ 附近抖动得厉害，绕数在不同网格下结果不一样。

原因：det(I+L)接近零时，行列式是多个因子的连乘，相对误差被放大，曲线在原点附近来回甩。smin小于1e-3*smax时必然出现。

解决：放弃行列式法，改用逐条特征值轨迹的相位法，并提升矩阵计算精度。如果smin在判定频段上始终很小，说明系统离不稳定边界很近，这时候讨论「稳不稳」本身意义有限，应该退回去改控制器。

5.2 判据层面：整数算对了，公式用错了

现象四：曲线绕数算成 1，系统开环也是稳定的，判据却说闭环不稳定，仿真却显示收敛。

原因：漏掉了镜像系数 2。只算了正频率半段的绕数，$Z=P+N$ 里的 $N$ 少了一半。

解决：在正频率段结果上乘 2，或者把判据写成 $Z=P+2N_+$。写代码时把这条放进函数名里，比如winding_positive_half，让调用方无法忽略这个语义。

现象五：开环存在右半平面极点，判据说稳定，实际系统发散。

原因：把 $P$ 当成了 0。对象本身不稳定时，$P>0$，此时曲线顺时针绕一圈对应的其实是有理函数零点极点抵消后的净值，不能孤立看 $N$。

解决：先用状态矩阵特征值或者分母求根算出 $P$，把它作为显式参数传进判据函数，不要写成默认值。

现象六：曲线在 $\omega=0$ 或者高频端出现「跳出去再回来」，绕数怎么算都不对。

原因：$L(s)$ 含积分环节或存在直通项，det(I+L)在端点发散，采样点在极限附近数值巨大，连线时直接跨过半个平面。

解决：把扫频区间从 $[0,\infty)$ 改成 $[\epsilon,\omega_{\max}]$，并单独解析处理小圆弧和大圆弧段；或者绕开频率响应，直接用状态空间模型求闭环极点做交叉验证。

6. 设计取舍：广义奈奎斯特、μ 分析与闭环特征值，谁该上场

四种手段经常被混着用，各自的保守性和代价差得很远。

小增益定理是充分条件里最容易验证的，把 $\sigma_{\max}(G K)$ 在整个频率范围压到 1 以下，一定稳，代价是要求极高，很多正常工作的系统都过不了。广义奈奎斯特是精确判据，但它的精确只针对标称模型：模型本身有误差时，曲线离 $-1$ 多近才算安全，它不回答。μ 分析回答的正是这个问题，代价是每个频点要做一次 D 尺度优化，而且上界和下界经常差出一截，得靠下界确认保守性到底有多大。

我的取舍习惯是分三层。设计阶段用广义奈奎斯特曲线看形状，找哪个模态、哪个频段在逼近 $-1$，这一步的反馈最快，图上直接看出裕度方向。定稿前用 μ 上界扫一遍，确认存在模型误差时裕度还在。最后用闭环状态矩阵的特征值做一次交叉验证，这一层最硬——它不依赖任何判据公式，纯粹是数值求根，如果它和曲线判定的结论矛盾，多半是曲线这一侧出了问题，去查网格、配对和 $P$。

交叉验证还有一个细节值得做：把闭环特征值代回 $\det(I+L(s))$ 算残差，残差应该接近零。这个检查能抓出零极点对消——如果某个闭环极点对应的 $\det$ 值明显不为零，说明分子分母里有公因子被消掉了，此时用 $\det$ 零点数闭环极点的前提就不成立，整条判据的可信度都要打折。

最后说个血泪教训。早些年我做一个三输入三输出的并联机构模型，用单回路方法一条一条调，三条回路的奈奎斯特曲线都离 $-1$ 远远的，Bode 图上裕度看着也漂亮，结果整机上电就低频抖。后来把三条特征轨迹一起画出来，才发现两条轨迹在 3 rad/s 附近几乎贴到了一起，行列式曲线正好从原点旁边擦过去，单看任何一条回路都发现不了。从那以后，我每次做 MIMO 的稳定性判定，都强制走一遍「特征值轨迹看形状、行列式曲线定绕数、闭环极点做交叉验证」这三步，一步都不跳。希望这些能帮到你，少走一遍我当年那个弯路——那份 matlab 实现里把扫频加密、轨迹配对和绕数校验串成了一条完整链路，正好可以对着源码把上面这几层逐个对照看。

本文还有配套的精品资源，点击获取