简介面向多输入多输出通信系统稳定性分析与教学验证的MATLAB代码包适合无线通信、信号处理方向的研究生、教师及工程技术人员使用。内容围绕广义奈奎斯特曲线绘制展开涉及开环传递函数矩阵建模、符号推导与模型降阶等环节帮助读者在复平面上判断系统稳定边界与相位特性。压缩包共7个文件均为.m脚本体积约5KB涵盖主绘图函数、示例脚本及配套的符号转换与降阶辅助函数结构紧凑便于逐行阅读与二次修改。目前已有3166人学习说明其在相关教学与科研场景中具备一定参考价值。读者可直接运行示例观察不同参数与信道条件下的曲线变化并对照源码理解各步骤的数学含义从而把理论公式、MATLAB编程与稳定性判据串联起来形成可复用的分析流程。1. 三个回路各自都稳合起来却发散了MIMO 稳定性判据缺的那块拼图三个各自稳定的单回路接进同一台多输入多输出设备之后整体振荡——这种场景在 matlab 里能复现在台架上也能复现。直觉上每条回路的开环奈奎斯特曲线都离 −1 点远得很闭环能差到哪去问题恰恰出在「单独」两个字上MIMO 系统的稳定性由整条特征轨迹共同决定断开其余回路去测其中一条测到的只是「其他回路还开着」的时候那个局部情况。多输入多输出本质是 m×m 的传递矩阵耦合项让每个输出同时受多个输入影响相位裕度和增益裕度就不再是标量而变成一组随频率旋转的数。想真正判断闭环稳不稳得把 m 条开环特征轨迹一并放到复平面上观察它们是否作为整体绕过了 −1 点。这套判据叫广义奈奎斯特判据围绕它画出来的曲线就是广义奈奎斯特曲线。这份 matlab 资源要解决的正是「怎么把这条曲线画对、并且从图上读出一个可信的整数结论」这件事。2. 广义奈奎斯特判据为什么盯住 det(IL) 而不是逐个回路2.1 闭环极点藏在 det(IL(s)) 的零点里单位负反馈下开环传递矩阵记作 $L(s)G(s)K(s)$$G$ 是被控对象$K$ 是控制器两者都是 $m\times m$ 的有理函数矩阵。闭环传递矩阵是 $(IL)^{-1}L$它的极点只可能来自两处$(IL)^{-1}$ 的奇点以及 $L$ 自身零极点对消后残留的部分。前者由 $\det(IL(s))0$ 决定后者在无对消假设下不贡献右半平面极点。于是闭环稳定性就被压缩成一个标量有理函数的零点位置问题$\det(IL(s))$ 在右半平面有没有零点。这一步转换是整个判据的地基也是它和单回路奈奎斯特的分水岭。单回路里判的是 $1L(s)$ 这条曲线绕不绕原点多回路里判的是 $\det(IL(s))$ 这条曲线绕不绕原点形式上一模一样含义从「一条回路」变成了「整个矩阵的联合行为」。把 $\det(IL(s))$ 写成 $N_c(s)/D_L(s)$$D_L$ 是 $L$ 各元素分母的最小公倍式$N_c$ 就是闭环特征多项式相差一个常数因子。开环不稳定极点个数 $P$ 等于 $D_L$ 在右半平面的根数闭环不稳定极点个数 $Z$ 等于 $N_c$ 在右半平面的根数。对 $\det(IL)$ 沿右半平面围线用幅角原理得到 $Z-PN$$N$ 是曲线顺时针绕原点的圈数。约定符号必须一开始就固定否则后面全乱。工程上习惯取「顺时针为正」因为顺时针绕圈对应右半平面零点多于极点。这个约定一旦写进代码注释就不要在中途改。2.2 乘积结构绕 −1 与绕原点是同一件事$\det(IL)$ 有一个非常好用的展开$\det(IL)\prod_{i1}^{m}(1\lambda_i(L))$$\lambda_i$ 是 $L$ 的特征值。每个因子 $1\lambda_i$ 在复平面上就是从 $-1$ 指向 $\lambda_i$ 的向量。所以 $\det(IL)$ 的相位等于各因子相位之和模等于各因子模之积。这个乘积结构直接给出两种等价的画法。第一种直接画 $\det(IL(j\omega))$ 这一条曲线数它绕原点的圈数判据是 $ZPN_{\det}$。第二种把 $m$ 条 $\lambda_i(L(j\omega))$ 轨迹画出来数它们各自绕 $-1$ 点的圈数之和判据是 $ZP\sum_i N_i$。两条路在数学上完全等价因为「某个 $1\lambda_i$ 转过原点」和「某个 $\lambda_i$ 转过 $-1$」是同一件事。等价不等于等效。行列式法的优点是只要一条曲线视觉上干净缺点是行列式是连乘任意一个因子接近零都会把整条曲线拽到原点附近数值相对误差被放大而且曲线绕圈时你无法判断是哪个模态干的。特征值法的优点是能定位到具体模态哪个回路、哪个方向在逼近 $-1$ 一目了然调控制器时这个信息很值钱缺点是特征值本身在重根附近不可微轨迹会出现真正的交叉连线逻辑要额外处理。对比项行列式法特征值轨迹法曲线条数1 条m 条数值稳定性近奇异时差单条轨迹较好排序依赖强模态可解释性无强重根处行为平滑可能交叉跳变适合场景快速判定稳定性调试、定位、做设计迭代我一般两个都画特征值轨迹用来看「谁在惹事」行列式曲线用来做最终判定两者绕数对不上就说明采样或数值出了问题。2.3 P 从哪来被控对象自身的右半平面极点$P$ 是 $\det(IL)$ 在右半平面的极点数在无零极点对消的前提下等于 $L$ 的右半平面极点数。$LGK$所以 $P$ 由被控对象和控制器共同决定。绝大多数人被这个量误导过以为单位反馈、控制器稳定$P$ 就是零。只有被控对象本身稳定、控制器也稳定时 $P0$ 才成立。一旦对象存在右半平面极点——某些机械结构、某些化工过程的线性化模型都有——$P$ 就非零此时哪怕曲线绕了 −1 一圈系统照样可能稳定因为 $ZPN$ 里 $P$ 把它抵消掉了。反过来$P0$ 而曲线顺时针绕了一圈那 $Z1$闭环必有一个右半平面极点。实际计算 $P$ 的稳妥做法是直接求 $L(s)$ 所有元素分母的根统计实部大于零的个数注意去掉重根。用状态空间模型更省事如果 $G$ 有状态矩阵 $A$$P$ 就是 $A$ 的特征值中实部为正的个数前提是 $G$ 的实现是最小实现。非最小实现里那些不可控不可观的模态会被算进 $P$导致判据整体偏移。提示$P$ 算错的后果比采样不够密严重得多因为它是一个整数的整体平移。先确认 $P$再谈曲线。3. 扫频与特征值配对奈奎斯特曲线能不能连对3.1 对数网格为什么会在谐振峰处骗你画频率响应习惯上取 $\omega$ 的对数等间距网格从 $10^{-2}$ 到 $10^{3}$ 取几百个点。这对 Bode 图够用对奈奎斯特曲线经常不够。原因是曲线在 $-1$ 附近的曲率由 $\sigma_{\min}(IL)$ 的极小值决定而在谐振峰附近特征值随频率的变化率可以大到相邻两个频点之间的相位增量超过 $\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 endC(i,j)是上一频点第 $i$ 条轨迹到当前频点第 $j$ 个特征值的欧氏距离greedy_assign做的是按行贪心的一一匹配。严格意义上这应该用匈牙利算法求全局最优但频率步长足够小时贪心解和最优解几乎总是一致代价小很多。Ci(used) inf那一行是防止两条轨迹抢同一个特征值丢掉它会在重根附近产生重复分配。这里有个容易绕进去的概念即使配对错了绕数在数学上也不该变。因为 $\lambda$ 集合的并集才是决定绕数的对象分支如何连接不影响整体是否包围 $-1$——前提是轨迹不穿过 $-1$。所以配对的作用主要是两件事让图能看让基于单条轨迹的算法比如逐条算相位不出错。真正做稳定性判定时用行列式法或对全部轨迹求和可以绕开配对问题。理解到这层就知道配对是工程手段不是数学必需。3.3 用最小奇异值给曲线可信度打分特征值在重根附近不可微数值上表现为两列特征向量接近共线eig的结果对扰动极其敏感。判断某个频点的特征值能不能信看 $\sigma_{\min}(IL(j\omega))$它是 $IL$ 到奇异矩阵的最小距离也近似是特征值 $\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 endsmin接近零的频点要单独标记出来这些点上det(IL)的相对误差会被放大到不可接受行列式曲线的位置不可信。经验阈值是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(dpi,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(k1); if (y1 0) ~ (y2 0) % 两端点在实轴两侧 xi x(k) - y1 * (x(k1)-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}(IL)$ 是否接近零而不是继续加压网格。还有一个频率端点的问题。当 $L(s)$ 含积分环节时$s0$ 是极点$\det(IL(j0))$ 发散正频率段不能从 0 开始扫得从 $10^{-4}$ 之类的下限起。当 $\det(IL(j\omega))$ 在 $\omega\to\infty$ 不趋于常数时高频端要做同样的截断。这两个截断不影响统计 RHP 包围次数的整数值——根据小圆弧或大圆弧的渐近方向包围次数要么不变要么只差一个可以解析确定的项。但很多画出来的曲线在极限处直接跳到无穷再回来肉眼看不出来这一跳对应多少次包围遇到这种情况宁可回到状态空间模型用精确的闭环极点核对也不要硬数图形。关于镜像正频率段的绕数要乘 2 才是完整的 RHP 围线绕数。因为实系数有理函数满足 $L(-j\omega)\overline{L(j\omega)}$负频率段是正频率段的共轭反向两者贡献相同相加就是加倍。这一步非常容易被漏掉漏掉之后 $ZPN$ 会把「不稳定」和「稳定」判反而且反得很有规律——你可能会发现模型明明发散判据却说稳定。5. 排查曲线画出来像那么回事结论却是反的5.1 图形层面曲线连得对眼睛才信得过现象一曲线在某一段出现横跨复平面的直线看起来像突然绕了一大圈。原因eig在相邻频点返回的顺序发生交换而代码里没做轨迹配对直接按特征值下标连线。重根附近这个问题几乎必然出现。解决在频率方向做最近邻配对见 3.2 的eigen_trajectories或者在绘图前用sort按实部排序——后者只是把突跳从随机位置挪到固定位置不能真正修好只能用于快速看一眼。正式判定必须做配对。现象二图上明显没绕过 −1判定结果却说绕了一圈。原因频率网格在谐振峰附近跨过了整段圆弧。曲线在两个采样点之间实际绕了 −1 一整圈又回来采样点上完全看不出来。解决用|1\lambda_i|的局部极小值定位关键频段在极小值附近用几何中点加密直到单步相位增量低于 $\pi/4$。更省事的办法是看smin曲线smin的谷底就是必须加密的地方。现象三曲线在 $-1$ 附近抖动得厉害绕数在不同网格下结果不一样。原因det(IL)接近零时行列式是多个因子的连乘相对误差被放大曲线在原点附近来回甩。smin小于1e-3*smax时必然出现。解决放弃行列式法改用逐条特征值轨迹的相位法并提升矩阵计算精度。如果smin在判定频段上始终很小说明系统离不稳定边界很近这时候讨论「稳不稳」本身意义有限应该退回去改控制器。5.2 判据层面整数算对了公式用错了现象四曲线绕数算成 1系统开环也是稳定的判据却说闭环不稳定仿真却显示收敛。原因漏掉了镜像系数 2。只算了正频率半段的绕数$ZPN$ 里的 $N$ 少了一半。解决在正频率段结果上乘 2或者把判据写成 $ZP2N_$。写代码时把这条放进函数名里比如winding_positive_half让调用方无法忽略这个语义。现象五开环存在右半平面极点判据说稳定实际系统发散。原因把 $P$ 当成了 0。对象本身不稳定时$P0$此时曲线顺时针绕一圈对应的其实是有理函数零点极点抵消后的净值不能孤立看 $N$。解决先用状态矩阵特征值或者分母求根算出 $P$把它作为显式参数传进判据函数不要写成默认值。现象六曲线在 $\omega0$ 或者高频端出现「跳出去再回来」绕数怎么算都不对。原因$L(s)$ 含积分环节或存在直通项det(IL)在端点发散采样点在极限附近数值巨大连线时直接跨过半个平面。解决把扫频区间从 $[0,\infty)$ 改成 $[\epsilon,\omega_{\max}]$并单独解析处理小圆弧和大圆弧段或者绕开频率响应直接用状态空间模型求闭环极点做交叉验证。6. 设计取舍广义奈奎斯特、μ 分析与闭环特征值谁该上场四种手段经常被混着用各自的保守性和代价差得很远。方法依据保守性计算代价能回答什么小增益定理$\sigma_{\max}(L)1$最高极低增益有界下的稳定性广义奈奎斯特$\det(IL)$ 或 $\lambda_i(L)$ 的绕数中对结构不确定性偏保守中需要加密扫频标称模型下闭环精确稳定与否结构奇异值 μ$\min_D \sigma_{\max}(DMD^{-1})$ 的上界低接近充要高反复做 D 尺度优化给定结构不确定性下的鲁棒稳定性闭环特征值$A_{cl}$ 的谱无只针对标称模型低标称模型极点位置小增益定理是充分条件里最容易验证的把 $\sigma_{\max}(G K)$ 在整个频率范围压到 1 以下一定稳代价是要求极高很多正常工作的系统都过不了。广义奈奎斯特是精确判据但它的精确只针对标称模型模型本身有误差时曲线离 $-1$ 多近才算安全它不回答。μ 分析回答的正是这个问题代价是每个频点要做一次 D 尺度优化而且上界和下界经常差出一截得靠下界确认保守性到底有多大。我的取舍习惯是分三层。设计阶段用广义奈奎斯特曲线看形状找哪个模态、哪个频段在逼近 $-1$这一步的反馈最快图上直接看出裕度方向。定稿前用 μ 上界扫一遍确认存在模型误差时裕度还在。最后用闭环状态矩阵的特征值做一次交叉验证这一层最硬——它不依赖任何判据公式纯粹是数值求根如果它和曲线判定的结论矛盾多半是曲线这一侧出了问题去查网格、配对和 $P$。交叉验证还有一个细节值得做把闭环特征值代回 $\det(IL(s))$ 算残差残差应该接近零。这个检查能抓出零极点对消——如果某个闭环极点对应的 $\det$ 值明显不为零说明分子分母里有公因子被消掉了此时用 $\det$ 零点数闭环极点的前提就不成立整条判据的可信度都要打折。最后说个血泪教训。早些年我做一个三输入三输出的并联机构模型用单回路方法一条一条调三条回路的奈奎斯特曲线都离 $-1$ 远远的Bode 图上裕度看着也漂亮结果整机上电就低频抖。后来把三条特征轨迹一起画出来才发现两条轨迹在 3 rad/s 附近几乎贴到了一起行列式曲线正好从原点旁边擦过去单看任何一条回路都发现不了。从那以后我每次做 MIMO 的稳定性判定都强制走一遍「特征值轨迹看形状、行列式曲线定绕数、闭环极点做交叉验证」这三步一步都不跳。希望这些能帮到你少走一遍我当年那个弯路——那份 matlab 实现里把扫频加密、轨迹配对和绕数校验串成了一条完整链路正好可以对着源码把上面这几层逐个对照看。本文还有配套的精品资源点击获取