☰
MIMO广义奈奎斯特曲线MATLAB绘制与绕数稳定判据
2026/11/2 13:35:46 网站建设 项目流程

简介:面向多输入多输出系统稳定性分析与频域可视化的 MATLAB 代码包,适合通信、控制与信号处理方向的学生和研究人员用来理解广义奈奎斯特曲线的计算与绘制。压缩包为 rar 格式,共 7 个文件,均为 .m 脚本,整体约 5KB,体积轻量却覆盖系统建模、传递函数符号转换、模型降阶、广义奈奎斯特曲线计算及示例验证等环节,便于从模型准备一路走到复平面绘图,形成从理论验证到结果观察的完整链路。已有 3166 人学习,说明其在 MIMO 频域分析入门与实验复现中具有一定参考价值。读者可借助这些脚本观察多输入多输出系统的稳定边界与相位特性,结合示例调整参数、模拟不同信道条件,进而加深对 MIMO 系统性能与稳定性判据的理解,也可作为课程实验或课题预研的辅助工具。

1. MIMO系统的奈奎斯特曲线:为什么不能把SISO的-1点判据直接搬过来

多输入多输出MIMO系统做频域稳定性校核,广义奈奎斯特曲线绘制绕不开一个动作:把开环传递矩阵L(jω)的特征值轨迹画出来,看它们对(-1,0)点的包围关系。很多刚接触matlab MIMO稳定性分析的人会直接把SISO的奈奎斯特曲线套到每个通道上,结果在耦合强的对象上翻车。这份资源是一套围绕MATLAB的脚本与示例,重点不是画一条漂亮曲线,而是把L(jω)特征值轨迹、围线闭合、绕数统计和稳定判据串起来。它适合做多变量频域校核、控制器增益试算和课程设计复现的人;如果你只画过单回路Nyquist,第一次用会碰到特征值排序、负频率镜像和积分器绕行这几个硬问题。

2. 建模与开环矩阵:把G(s)K(s)整理成可扫频的L(s)

2.1 先定开环回路:L(s)=G(s)K(s)而不是G(s)

广义奈奎斯特判据看的是闭环回差矩阵I+L(s),其中L(s)是开环传递矩阵。MIMO系统里,G(s)是对象,K(s)是控制器,常见做法是先算L=GK,再对L(jω)做特征值分解。这里最容易出的错是把G(s)直接拿去画图,然后在图上找不到-1点对应的包围关系。SISO里开环是g(s)k(s),MIMO里矩阵乘法不交换,GK和K*G的特征值一般不同,回差矩阵也不一样。脚本里要固定一种回路定义,后面P统计、绕数统计和时域仿真都用同一个L,否则判据结果会和仿真对不上。

我一般先确认三件事:对象是严格真还是双真、控制器有没有积分、L(s)在s=0处有没有极点。严格真意味着L(∞)=0,高频段特征值轨迹会往原点收;双真意味着L(∞)=D,高频段轨迹会停在D的特征值上,围线闭合时无穷远大半圆映射成一个点而不是原点。控制器带积分时,L(s)在s=0有极点,freqresp(L,0)可能返回Inf,扫频不能从0直接开始,要在右半平面做小半圆绕行。这些不是绘图细节,是判据成立的前提。

资源里的示例通常会给一个2x2对象和一个对角增益矩阵,方便先跑通。拿到脚本后不要急着换自己的模型,先用示例确认坐标比例、-1点位置和绕数统计是否正常,再替换G和K。这个顺序能省掉很多“图看起来对、判据结果不对”的排查时间。

2.2 tf与ss的取舍:freqresp返回的维度要先验

MATLAB里搭L(s)可以用tf、ss或zpk。做频域扫频时,tf写起来直观,ss在状态维数高时数值更稳,zpk适合零极点清楚的模型。对广义奈奎斯特绘制,关键不是用哪种表示,而是freqresp返回的数组维度要搞清楚。MIMO的freqresp(sys,w)返回的是[ny,nu,length(w)]的三维数组,ny是输出数,nu是输入数,第三维是频率点。后面eig要逐频点取ny×nu矩阵,不能把三维数组直接丢给eig。

下面这段代码先构造一个2x2开环L,再验证freqresp的维度。示例对象用一阶惯性加耦合项,控制器用对角增益,目的是让特征值轨迹有交叉但不至于一上来就太复杂。

% 2x2 MIMO开环模型:G(s)K(s) s = tf('s'); G = [1/(s+1), 0.6/(s+2); 0.4/(s+3), 1/(s+1.5)]; K = diag([1.2, 0.8]); % 对角控制器,先保证回路定义清楚 L = G * K; % 先选一段正频率,验证freqresp维度 w_test = logspace(-2, 2, 5); Lw_test = freqresp(L, w_test); size(Lw_test) % 期望:2 2 5

逻辑说明:L=G*K是开环回路,后续所有特征值、绕数和稳定判据都基于这个L。freqresp不直接返回传递矩阵的符号表达式,而是返回每个频点的复数值矩阵,所以第三维长度必须和频率向量长度一致。参数说明:w_test只是维度检查用的粗网格,真正绘图时频率点要加密到几千点;K先取对角矩阵,是为了让L的特征值和通道增益有直观对应,等流程跑通再换非对角控制器。

如果size返回的不是[2,2,5],先检查G和K是不是都定义成了LTI对象,以及有没有把某个通道误写成标量。另一个常见问题是把freqresp的返回值直接当二维矩阵用,比如Lw_test(:,:,1)才是第一个频点的2x2矩阵。这个维度错误在绘图阶段往往表现为“eig输入不是方阵”或者“只画出一条线”,排查时先打印size比盯着图看有效。

2.3 频率网格与围线形状:logspace之外的三个参数

频率网格不是越密越好,但广义奈奎斯特曲线对-1点附近的穿越很敏感,网格太稀会把包围数数错。常见做法是用logspace覆盖低频到高频,低频下限根据最慢模态或积分器绕行半径定,高频上限根据带宽、增益穿越频率和L(∞)是否为零定。如果L(s)严格真,高频段轨迹会收向原点,上限取到比对象带宽高两到三个数量级通常够用;如果L(s)双真,高频段停在D的特征值上,上限不需要无限拉高,但要把D的特征值看清楚。

第二个参数是频率点密度。正频率段用6000到10000点不算多,尤其是特征值轨迹在-1点附近拐弯时,稀疏网格会把绕数算成0。第三个参数是负频率镜像。实系数系统的负频率响应是正频率的共轭,画完整围线时不能只画正频率,否则曲线不闭合,绕数统计没有意义。下面给一个常用的频率网格构造方式,正频率从一个小正数开始,负频率用共轭镜像补全。

% 频率网格:正频率 + 负频率镜像 w_lo = 1e-4; % 低频下限,有积分器时不要取0 w_hi = 1e4; % 高频上限,按对象带宽调整 Nw = 8000; % 正频率点数 w_pos = logspace(log10(w_lo), log10(w_hi), Nw); % 负频率按共轭镜像补全:-w_hi ... -w_lo w_full = [w_pos, -fliplr(w_pos)]; % 如果L没有虚轴极点,可以在两端补0附近点让曲线更闭合 % 有积分器时不要直接补0,改用右半平面小半圆绕行

逻辑说明:w_pos负责正频率扫频,w_full把负频率按从高到低镜像接在后面,形成一条从低频到高频、再从负高频回到负低频的频率序列。参数说明:w_lo不能取0,因为logspace(0,...)会出问题,而且有积分器时freqresp在0处可能无穷大;w_hi按你的执行机构带宽和对象高频衰减定,如果高频段轨迹已经缩成一点,再往上加只会增加计算量;Nw影响绕数统计的稳定性,第一次跑可以用8000,后续用加密网格做交叉检查。

提示:频率网格和围线形状要一起看。只改网格不改围线,负频率漏画时曲线仍然不闭合;只补负频率不加密网格,-1点附近的穿越仍可能数错。

3. MATLAB扫频绘制广义奈奎斯特曲线:特征值轨迹与完整围线

3.1 特征值轨迹计算:eig逐频点与最近邻排序

MIMO广义奈奎斯特曲线的核心是L(jω)的特征值轨迹。每个频点对L(jω)做eig,得到n个特征值,n是方阵维数。问题在于eig每次返回的特征值顺序不保证连续,频率稍微一变,两个特征值可能互换位置,画出来的线就会突然跳。这个现象在特征值接近或交叉时特别明显,图上看起来像“断线”或者“多出一条莫名其妙的轨迹”,实际是排序问题。

常见做法是逐频点做最近邻匹配:当前频点的每个特征值,去上一频点的特征值里找距离最近的那个,按这个对应关系排。这个方法不完美,特征值真正交叉时最近邻也会选错,但配合加密扫频和后面的绕数交叉检查,足够把大多数工程问题跑通。更稳的做法是用子空间跟踪或匹配算法,但脚本里先用最近邻能把流程搭起来。

function lam = eig_track(Lw) % 输入Lw: ny x nu x N,且ny=nu % 输出lam: n x N,按相邻频点最近邻排序的特征值轨迹 n = size(Lw,1); N = size(Lw,3); lam = zeros(n,N); lam(:,1) = eig(Lw(:,:,1)); for k = 2:N ek = eig(Lw(:,:,k)); prev = lam(:,k-1); used = false(n,1); order = zeros(n,1); for j = 1:n d = abs(ek(j) - prev); d(used) = inf; % 已经被匹配的上一频点特征值不再参与 [~, idx] = min(d); order(j) = idx; used(idx) = true; end lam(:,k) = ek(order); end end

逻辑说明:先算第一个频点的特征值作为基准,后面每个频点都拿当前特征值和上一频点已排序的特征值做距离匹配,距离最近的一对认为属于同一条轨迹。参数说明:Lw必须是方阵频响,ny=nu;如果系统不是方阵,广义奈奎斯特特征值轨迹的定义要先做方阵化处理,不能直接套这个函数;距离用复平面欧氏距离abs(ek(j)-prev),没有做归一化,增益跨度大时可以先对L做缩放,后面避坑章节会讲。

这个函数跑完后,lam的第i行就是第i条特征值轨迹。绘图时对每一行画real和imag。要注意:最近邻排序只在相邻频点之间成立,频率网格太稀时上一频点和当前频点之间可能已经转过大半圈,最近邻会匹配到错误的轨迹。第一次跑不要用太稀的网格,看到轨迹跳线先加密频率,再怀疑排序算法。

3.2 补全负频率与无穷远圆弧:闭合曲线才有绕数

广义奈奎斯特判据数的是闭合围线的包围数,只画正频率轨迹只能看趋势,不能直接数绕数。实系数系统的负频率响应是正频率的共轭,所以负频率轨迹可以用conj(fliplr(lam))补出来。频率序列也要对应成[w_pos, -fliplr(w_pos)],这样曲线从低频走到高频,再从负高频走回负低频,形成一条接近闭合的路径。

如果L(s)严格真,高频段L(jω)→0,正高频和负高频在原点附近接上,无穷远大半圆映射到原点,曲线自然闭合。如果L(s)双真,L(∞)=D,高频段停在D的特征值上,正负频率在D的特征值处相接,无穷远大半圆映射成D的特征值点,不是一个圆弧。积分器对象另说,s=0附近要做右半平面小半圆绕行,不能直接用0频点。

下面代码在上一节的eig_track基础上补全负频率并绘图。为了让-1点看得清楚,坐标轴用等比例,-1点用醒目标记。

% 正频率扫频 w_pos = logspace(-4, 4, 8000); Lw_pos = freqresp(L, w_pos); lam_pos = eig_track(Lw_pos); % n x N % 负频率共轭镜像:注意顺序要反过来 lam_neg = conj(fliplr(lam_pos)); lam_full = [lam_pos, lam_neg]; % 绘图 figure; hold on; grid on; for i = 1:size(lam_full,1) plot(real(lam_full(i,:)), imag(lam_full(i,:)), ... 'LineWidth', 1.2); end plot(-1, 0, 'kx', 'MarkerSize', 12, 'LineWidth', 2); % -1点 axis equal; xlabel('Re(\lambda(L(j\omega)))'); ylabel('Im(\lambda(L(j\omega)))'); title('MIMO广义奈奎斯特曲线:L(j\omega)特征值轨迹');

逻辑说明:lam_full把正频率轨迹和负频率共轭镜像拼在一起,形成完整围线的映射。plot逐条画特征值轨迹,-1点用黑色叉号标出,因为MIMO判据看的是特征值轨迹对-1点的包围,不是对原点的包围。参数说明:w_pos的范围按上一节的w_lo和w_hi调整;lam_neg用conj(fliplr(...)),fliplr保证负频率从高频往低频走,conj保证实系数系统的共轭对称;axis equal必须加,否则坐标轴比例失真,包围关系看起来会变形。

如果曲线在低频端没有闭合,先看w_lo是不是太大,或者L(s)在低频有极点需要绕行。如果曲线在高频端没有闭合,看L(s)是不是严格真,严格真时高频轨迹应该收向原点;如果停在某个非零点,检查L(s)是不是双真,双真时高频停在D的特征值上是正常的。不要为了“看起来闭合”强行把首尾用直线连起来,那条直线在s平面没有对应的围线映射,数出来的绕数没有意义。

3.3 绘图细节:-1点、箭头方向与坐标等比例

图画出来只是第一步,读图时最容易错的是方向。Nyquist围线在s平面是顺时针包围右半平面,映射到L平面后,特征值轨迹的方向决定了包围数是正还是负。脚本里可以在几个频点加箭头,标出ω增加的方向。MATLAB没有直接的箭头函数,常见做法是用quiver在相邻点之间画小箭头,或者用annotation,但annotation坐标是归一化的,换图就乱。我一般用quiver在轨迹上隔一段画一个箭头,简单直接。

另一个细节是-1点的位置。MIMO判据看的是λ_i(jω)对-1点的包围,不是det(I+L)对原点的包围。两者等价,但图上标的点不同。如果你在图上标了原点又数-1点的包围,肯定对不上。绘图时把-1点标成醒目的叉号,原点可以标一个淡一点的点作为参考,但判据只认-1点。

坐标范围也要注意。特征值轨迹可能在低频段有很大的实部或虚部,自动坐标轴会把-1点附近压得很小,包围关系看不清。常见做法是先画全图,再用xlim/ylim缩放到-1点附近做局部检查。局部图不能替代全图数绕数,但能看清轨迹是从-1点上方还是下方穿过。下面给一个局部放大和箭头标注的片段。

% 在完整轨迹上每隔一段加方向箭头 step = round(size(lam_full,2)/30); for i = 1:size(lam_full,1) for k = 1:step:size(lam_full,2)-1 x0 = real(lam_full(i,k)); y0 = imag(lam_full(i,k)); dx = real(lam_full(i,k+1)) - x0; dy = imag(lam_full(i,k+1)) - y0; quiver(x0, y0, dx, dy, 0, 'Color', [0.4 0.4 0.4]); end end % 局部放大到-1点附近,检查穿越方向 xlim([-2.5 0.5]); ylim([-1.5 1.5]);

逻辑说明:quiver在相邻频点之间画箭头,方向就是ω增加的方向;step控制箭头密度,太密会糊成一片,太稀看不出方向。xlim/ylim把视图缩到-1点附近,方便检查轨迹从上方还是下方穿过。参数说明:箭头颜色用灰色,避免盖住轨迹颜色;局部范围按你的模型调整,如果-1点附近轨迹很密,可以再放大到[-1.5,-0.5]×[-0.5,0.5];局部图看完要回到全图确认闭合,不能只看局部就数绕数。

4. 稳定判据怎么读图:净包围数、P统计与奇异值旁证

4.1 净包围数不是数穿越:用unwrap相位算绕数

广义奈奎斯特判据的结论可以写成:设开环右半平面极点数为P,L(jω)特征值轨迹对-1点的净逆时针包围数为N,则闭环右半平面极点数Z=P-N,稳定要求Z=0,即N=P。这里的N是净包围数,不是穿越次数。轨迹从-1点上方穿过去一次、再从下方穿回来一次,净包围数可能是0,但穿越次数是2。只数穿越会把条件稳定系统判错。

数值上算绕数,常用方法是看1+λ_i(jω)的相位沿闭合围线的总变化。对每条特征值轨迹,构造z_i=λ_i+1,计算angle(z_i)并unwrap,闭合围线的总相位变化除以2π就是绕数。方向符号和围线方向有关,第一次跑用已知稳定和不稳定的例子标定一下正负号,不要硬记。下面函数输入一条闭合的特征值轨迹,输出绕数。

function N = winding_number_track(lam_track) % lam_track: 一条闭合的特征值轨迹,首尾应接近同一点 % 返回对-1点的净包围数,逆时针为正(按你的围线方向标定) z = lam_track(:) + 1; % 平移到原点:λ对-1的包围 -> λ+1对原点的包围 theta = unwrap(angle(z)); N = (theta(end) - theta(1)) / (2*pi); end

逻辑说明:λ_i对-1点的包围等价于λ_i+1对原点的包围,所以先把轨迹平移1。angle取每个点的相角,unwrap展开跳变,闭合围线的总相角变化除以2π就是绕数。参数说明:lam_track必须是闭合序列,首尾要接上;如果首尾不严格相等,可以复制第一个点补到末尾,但补点只影响闭合,不影响轨迹形状;方向符号用你的围线方向标定,SISO里找个开环稳定且闭环稳定的例子,绕数应为0,再找个闭环不稳定的例子看符号。

把所有特征值轨迹的绕数加起来,得到总N。总N和det(I+L)的绕数一致,因为det(I+L)=∏(1+λ_i)。如果单条轨迹的绕数不好数,可以用det(I+L)的相位变化做交叉检查,最后一章会给代码。这里先记住:单条轨迹绕数和总绕数都要能对上,对不上时先查特征值排序,再查负频率有没有漏。

4.2 开环右半平面极点P:pole统计与积分器绕行

P是开环L(s)在右半平面的极点数。MIMO系统的极点不是每个通道极点简单相加,而是传递矩阵的最小实现极点。用pole(L)可以返回L的极点,统计实部大于0的个数就是P。注意数值误差:实部可能是1e-12这种接近0的数,统计时要设一个容差,比如real(p)>1e-8。积分器在虚轴上,实部为0,不算右半平面极点,但围线经过虚轴时要绕行,否则freqresp在极点处会发散。

积分器绕行的常见做法是在s=0附近取右半平面小半圆,s=ε e^{jθ},θ从-π/2到π/2,ε取一个很小的正数。这样围线不穿过极点,而是从右侧绕过。对应到频率响应,不能只用jω,要在小半圆上取复频率。MATLAB里freqresp支持复频率点,可以把s=ε e^{jθ}代入。下面片段演示统计P和构造小半圆频率点。

% 统计开环右半平面极点 p = pole(L); tol = 1e-8; P = sum(real(p) > tol); % 如果L在虚轴上有极点(如积分器),构造右半平面小半圆绕行 epsilon = 1e-4; theta = linspace(-pi/2, pi/2, 200); s_arc = epsilon * exp(1j*theta); % 右半平面小半圆 L_arc = freqresp(L, s_arc); % 注意:freqresp支持复频率

逻辑说明:pole(L)返回开环极点,real(p)>tol统计右半平面极点;小半圆用s=ε e^{jθ}从下方绕到上方,θ从-π/2到π/2,保证围线从右侧绕过虚轴极点。参数说明:tol不要设得太大,否则会把靠近虚轴的稳定极点算成右半平面;epsilon要小于最靠近虚轴的右半平面零点/极点的距离,太大改变围线形状,太小数值误差大,常见取1e-5到1e-3;theta点数200通常够,小半圆映射后变化剧烈时可以加到1000。

注意:积分器绕行不是可选项。L(s)有1/s时,freqresp(L,0)会返回Inf,直接扫到0会让整条轨迹出现NaN,后面的绕数统计全废。要么从w_lo>0开始并单独补小半圆,要么在围线里显式绕行。

4.3 奇异值曲线旁证:特征值轨迹什么时候会骗人

特征值轨迹是广义奈奎斯特判据的直接工具,但它对特征值排序和数值条件很敏感。MIMO系统在特征值接近时,轨迹会交叉,排序算法可能把两条轨迹接错,画出来的包围数看着像对的,实际是两条轨迹各转了一半。另一种情况是L(jω)的条件数很大,特征值对扰动很敏感,图上-1点附近稍微改一点增益,轨迹就跳过去,判据结果跟着变。

奇异值曲线可以做旁证。σ_min(I+L)接近0表示回差矩阵接近奇异,闭环接近不稳定;σ_min(I+L)在整个频率范围内的最小值可以当作稳定裕度的粗略指标。注意奇异值和特征值不是一回事:特征值轨迹看的是包围关系,奇异值看的是增益放大,两者不能互相替代,但可以互相印证。如果你在特征值轨迹上看到-1点附近有穿越,同时σ_min(I+L)在该频率附近明显下凹,那这个穿越大概率是真实的;如果特征值轨迹看着穿过了-1点,但σ_min(I+L)一直很大,先怀疑特征值排序或负频率漏画。

下面表格把三个量的用途分开,避免混用。

量看什么适合回答的问题容易误用的地方
λ_i(L(jω))对-1点的包围闭环稳定与否、绕数是否满足判据特征值排序跳线会被当成真实穿越
det(I+L(jω))对原点的包围总绕数交叉检查行列式接近0时相位跳变,需要unwrap
σ_min(I+L(jω))最小奇异值下凹稳定裕度旁证、接近奇异程度不能替代绕数判据,只做辅助

实际排查时,我一般先看特征值轨迹全图,再用det(I+L)绕数复核,最后看σ_min(I+L)在穿越频率附近有没有下凹。三步都对上,才敢说判据结果可信。只看一张图就下结论,在耦合强的对象上很容易翻车。

5. 避坑与排查:围线漏点、特征值排序、增益缩放的五条血泪记录

5.1 现象:闭合后包围数和det(I+L)对不上

现象:特征值轨迹图看起来闭合了,单条轨迹数出来的绕数是1,但det(I+L)的相位绕数是0,两个结果打架。原因通常是负频率镜像的顺序错了,或者正负频率在首尾没有接上。实系数系统的负频率响应是正频率的共轭,但顺序要从高频往低频走,才能和正频率的高频端接上。如果直接conj(lam_pos)不fliplr,负频率段的方向反了,闭合围线在-1点附近多绕了一圈,单条轨迹的绕数就会多1。

解决:检查lam_full的拼接顺序,必须是[lam_pos, conj(fliplr(lam_pos))]。频率序列也要对应[w_pos, -fliplr(w_pos)]。如果L(s)严格真,高频端应该收向原点,正高频和负高频在原点附近接上;如果L(s)双真,高频端停在D的特征值上,正负频率在D的特征值处相接。接不上时不要用直线强行连,先查L(∞)和频率上限。

另一个常见原因是低频端没有闭合。w_lo取1e-4时,正频率最低点和负频率最低点之间差一个0附近的小弧。如果L(s)没有虚轴极点,可以补一个0频点;如果有积分器,必须用小半圆绕行。漏掉这个小弧,绕数可能差半圈,表现出来就是单条轨迹和det(I+L)对不上。

5.2 现象:特征值轨迹突然跳线

现象:图上某条特征值轨迹在某个频率突然跳到另一条轨迹上,形成一条横跨的直线,或者两条轨迹交换位置后继续走。原因有两个:eig返回顺序不连续,或者频率网格太稀,相邻频点之间特征值已经转过大半圈,最近邻匹配到了错误的轨迹。前者加密网格能缓解,后者需要检查频率范围是不是覆盖了特征值的快速变化段。

解决:先用更密的频率网格重跑,比如从8000点加到20000点,看跳线是否消失。如果跳线还在,打印跳线频点附近的特征值,看是真的交叉还是排序错。真正交叉时,两条轨迹在复平面上相交,最近邻无法区分,可以改用子空间跟踪或匹配算法;排序错时,当前频点的特征值和上一频点的特征值距离并不近,最近邻却因为used标记选错了。脚本里可以加一个距离阈值,超过阈值时给出警告,提示该频点需要人工检查。

血泪经验:不要看到跳线就手动把两条轨迹交换回来。手动交换只对当前图有效,换一个频率网格又乱了。先把排序算法和网格密度调稳,再谈绘图美观。

5.3 现象:高频段一团乱麻、-1点附近全是线

现象:画出来的图高频段轨迹密密麻麻,-1点附近看不清哪条线穿过,绕数根本没法数。原因通常是频率范围太宽,把大量高频衰减段画进了同一张图,或者增益缩放不合适,-1点附近被压缩成一个小区域。MIMO系统的特征值轨迹在高频段可能收向原点,也可能停在D的特征值上,如果D的特征值离-1点很近,高频段会和判据区域混在一起。

解决:先分段画。低频段看积分器绕行和慢模态,中频段看穿越,高频段单独看L(∞)落在哪里。绘图时用xlim/ylim局部放大-1点附近,但绕数统计仍然用全图数据。增益缩放也要注意:如果L的增益很大,特征值轨迹整体远离-1点,图上看着稳定,但减小增益后可能穿过-1点;如果L的增益很小,轨迹挤在原点附近,-1点看起来很远,实际可能已经包围。常见做法是保持L的物理增益不变,用局部放大检查,不要为了图好看把L乘以一个系数。

5.4 现象:带积分器对象在低频报错或飞出画面

现象:L(s)里有1/s,freqresp(L,0)返回Inf或NaN,绘图时整条轨迹消失;或者低频段轨迹飞到很大的实部/虚部,坐标轴被拉爆,-1点附近完全看不见。原因是围线直接经过了虚轴极点,频率响应在极点处发散。MIMO系统的积分器可能出现在某个通道,也可能以矩阵形式出现,不能只看单个通道。

解决:在s=0附近做右半平面小半圆绕行,s=ε e^{jθ},θ从-π/2到π/2,把这段复频率响应单独算出来,接到正频率低频端和负频率低频端之间。ε取1e-5到1e-3,按最靠近虚轴的极点距离调。绘图时可以把小半圆段用不同颜色标出来,提醒自己这段是绕行段,不是真实频率响应。坐标轴不要用自动,先算全图范围,再手动设一个能看清-1点的范围,低频发散段可以截断显示,但绕数统计要用完整数据。

5.5 现象:判据说稳定,时域仿真却发散

现象:特征值轨迹算出来的绕数满足N=P,判据说闭环稳定,但时域仿真发散。原因可能是P统计算错、L的定义和仿真不一致、或者特征值排序把包围数算错了。P统计最容易在数值容差上出错:pole(L)返回的极点实部可能是1e-10,被算成右半平面,P多了1,判据要求的N也多了1,结果把不稳定判成稳定。另一个原因是仿真用的L和绘图用的L不是同一个回路,比如仿真里控制器有饱和或延迟,频域模型没包含。

解决:先核对P,用real(p)>1e-8统计,并把靠近虚轴的极点列出来人工看。再核对L,确保绘图和仿真用的是同一个G*K,控制器增益、符号、通道顺序都一致。最后用det(I+L)绕数交叉检查特征值轨迹的绕数,两者对不上时优先信det(I+L),因为它不依赖特征值排序。如果三者都对上但仿真仍发散,检查频域模型是否漏了延迟、采样或非线性,广义奈奎斯特判据只对线性时不变模型成立。

6. 进阶验证:用det(I+L)绕数与加密扫频交叉检查

6.1 用det(I+L)相位绕数做独立复核

特征值轨迹的绕数依赖排序,det(I+L)的绕数不依赖排序,适合做独立复核。做法是沿完整围线算det(I+L(jω)),取angle并unwrap,总相位变化除以2π就是绕数。这个绕数和所有特征值轨迹绕数之和应该一致。下面代码沿用前面的w_full和L,注意det(I+L)在行列式接近0时相位会跳,unwrap之后要检查有没有异常大的跳变。

% 用det(I+L)沿完整围线做绕数复核 F = zeros(1, length(w_full)); for k = 1:length(w_full) Lk = freqresp(L, w_full(k)); F(k) = det(eye(size(Lk,1)) + Lk); end theta = unwrap(angle(F)); N_det = (theta(end) - theta(1)) / (2*pi); % 如果首尾没有严格闭合,把第一个点补到末尾再算一次 F_closed = [F, F(1)]; theta_closed = unwrap(angle(F_closed)); N_det_closed = (theta_closed(end) - theta_closed(1)) / (2*pi);

逻辑说明:F是det(I+L)沿围线的复数值,angle取相角,unwrap展开跳变,总相位变化除以2π得到绕数。F_closed把首尾补上,避免围线不严格闭合导致的半圈误差。参数说明:w_full要和特征值轨迹用的频率序列完全一致;eye(size(Lk,1))保证回差矩阵维数正确;N_det的符号和围线方向有关,用已知稳定例子标定。

6.2 三个验证样例与失败信号

我一般用三个样例交叉检查:一个开环稳定且闭环稳定的2x2对象,绕数应为0;一个开环稳定但增益过大导致闭环不稳定的对象,绕数应为-1或+1(按方向标定);一个带积分器的对象,检查小半圆绕行后绕数是否连续。三个样例都通过,再换自己的模型。失败信号和对策可以列成表。

验证项期望失败时先看什么
特征值轨迹总绕数与det(I+L)绕数一致负频率拼接顺序、特征值排序
已知稳定样例绕数满足N=PP统计容差、L回路定义
已知不稳定样例绕数比稳定时多1-1点局部放大、频率网格密度
积分器样例小半圆段连续、无NaNε大小、绕行方向

加密扫频是最后一道后悔药。把Nw从8000加到20000,如果绕数变了,说明原网格漏掉了-1点附近的快速变化;如果绕数不变,结果基本可信。我现在的习惯是每次换模型先跑已知稳定和不稳定两个样例,再跑自己的对象,最后用det(I+L)复核一遍。从那以后我每次画广义奈奎斯特曲线,都强制走一遍“特征值排序检查→负频率闭合→det绕数复核”的流程,不再只看图下结论。希望帮到你。

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询