企业官网建设流程全解析

简介：面向无线通信与信号处理方向的研究人员、高年级本科生，提供一套基于MATLAB环境的多输入多输出系统广义奈奎斯特曲线绘制方案。借助 nyqmimo 函数与配套脚本，可针对开环传递函数矩阵快速生成广义奈奎斯特图，用于分析多输入多输出系统的稳定性边界、频率响应与相位特性，解决多变量耦合系统难以手工绘图和验证的难题。压缩包内共7个 M 文件，整体仅5KB，涵盖广义奈奎斯特曲线主绘制函数、符号传递函数转换、模型降阶以及可运行示例脚本等，代码量小、层次清楚，便于按需调用、修改或嵌入自身仿真流程。代码已经实测可用，改变系统参数或模拟不同信道条件时，可直观观察复平面曲线变化，进而深入理解广义奈奎斯特判据、系统稳定性与频域特性；同时也可作为课程设计、毕业设计或科研探索的辅助工具。目前已有3162人学习下载，是一份轻量、实用，对无线通信、信号处理教学与项目实践均有参考价值的学习与验证资料。

1. 多输入多输出MIMO系统的稳定性，为什么要画广义奈奎斯特曲线

多输入多输出（MIMO）系统的频域稳定性校验，最直接的工具就是广义奈奎斯特（Generalized Nyquist）曲线：把开环返回比矩阵 L(s)=G(s)K(s) 的特征值随频率 ω 在复平面上扫出的轨迹画出来，看它绕临界点 (−1, j0) 的总圈数。某次我在调一套双输入双输出温度回路的时候，两个通道单独看奈奎斯特判据都合格，合上耦合通道之后闭环就地振荡，最后就是靠特征轨迹定位到耦合支路的相位短板才救回来。

标量系统里只要画一条 G(s)K(s) 频率特性曲线就能判断稳定性，MIMO 回路里 G(s)K(s) 是一个矩阵，矩阵没有“相角”这个概念，稳定性信息藏在矩阵的特征值里。把返回比矩阵 L(s) 逐频点做特征值分解，每个特征值 λ_i(jω) 在复平面上形成的轨迹就是特征轨迹（characteristic loci），这些轨迹共同构成广义奈奎斯特曲线。

这篇笔记面向两类人：一类是刚接触多变量控制、想做频域校验的工程师，另一类是已经有 SISO 频域经验、被 MIMO 耦合绕晕的熟手。你会看到从数学构造到排序算法、绘图代码、判读方法、避坑清单，以及最后怎么把曲线结论送回时域仿真验证的完整路径。

2. 特征轨迹的数学构造：从闭环行列式到逐点特征值分解

2.1 为什么SISO奈奎斯特判据不能直接搬到MIMO

SISO 回路里，闭环特征方程是 1 + G(s)K(s) = 0，奈奎斯特曲线绕临界点 (−1, j0) 的净圈数直接决定右半平面闭环极点数。MIMO 回路里，闭环特征方程变成 det(I + G(s)K(s)) = 0，这里的 G(s)K(s) 是 m×m 矩阵，不能直接“画一条曲线”。

但矩阵的特征值分解能帮上忙。对任意方阵 A，特征值组成的对角矩阵为 Λ，满足 det(I + A) = det(I + Λ) = ∏(1 + λ_i)。把 A 换成 L(s) = G(s)K(s)，就有：

det(I + L(s)) = ∏ (1 + λ_i(s))

闭环极点就是 det(I + L(s)) = 0 的根，也就是任意一条特征轨迹满足 λ_i(s) = −1 的位置。所以整个 MIMO 稳定性问题被拆成了若干条标量轨迹的行为，每一条轨迹看它是否绕 (−1, j0)，统计总圈数即可。这里的 λ_i(s) 不是某个物理通道的开环增益，而是整个回路矩阵的特征增益，它天然包含了通道间耦合的影响。

反直觉的地方在于：单通道都稳定不代表矩阵稳定，因为特征增益的轨迹是耦合后的整体表现。两个通道各自的 G11 和 G22 相位裕量很足，交叉耦合 G12、G21 会把特征增益拉向临界点，这种情况只有画特征轨迹才看得到。

2.2 广义奈奎斯特判据的表达与符号约定

广义奈奎斯特判据（GNC）的工程表述并不复杂，但要先把符号约定写死，否则后面判读很容易出偏差。

判据公式是 Z = P + N。绝大多数实际工程里 P = 0，这时只要看 N 是否等于 0；如果开环回路里有右半平面极点，稳定条件变成 N = −P，也就是需要顺时针绕 P 圈。注意工程上常犯的错是把 P 只算被控对象 G(s) 的极点，忘了控制器 K(s) 也可能带来右半平面极点，这个坑在第 5 章会专门展开。

还有一个必须交代的约定：实有理传递矩阵满足 L(−jω) = conj(L(jω))，特征值也对应共轭。也就是说负频率段的特征轨迹不需要单独计算，直接把正频率段轨迹关于实轴做镜像就能补全。绘制时补全负频支，轨迹才能形成闭合曲线，方便用角度 unwrap 统计净圈数。

2.3 在频域逐点计算L(jω)的特征值

广义奈奎斯特曲线的数据来源一般有两种：解析传递函数矩阵和实测频率响应数据。先看解析模型怎么做，下面是一个典型的 2×2 被控对象，控制器用对角增益，按频点逐个求特征值。

import numpy as np def loop_TF2x2(s): # 2×2 开环对象： # 对角通道主增益 8/((s+1)(s+2)) # 耦合通道增益 1/(s+1) d = (s + 1) * (s + 2) return np.array([ [8.0 / d, 1.0 / (s + 1)], [1.0 / (s + 1), 8.0 / d], ]) def compute_eigen_curves(w_rad, kp): # 控制器为对角增益 kp*I，逐频点求返回比矩阵的特征值 vals = [] for w in w_rad: s = 1j * w L = loop_TF2x2(s) * kp # L(s) = G(s) * K(s) vals.append(np.linalg.eigvals(L)) return np.array(vals)

代码逻辑很直接：把 s = jω 代入传递函数矩阵，乘以控制器增益，然后用 numpy.linalg.eigvals 求特征值。关键点是 L(s) 的乘法顺序，这里是 G(s)K(s)，控制器在对象后面；如果是输出反馈、控制器在对象前面，L(s) = K(s)G(s)，特征值集合不变，但书写要一致。

频点向量 w_rad 的单位是 rad/s。如果系统穿越频率在 1 到 100 rad/s 之间，建议起点 0.001、终点 1000，跨越 6 个十倍频程，点数至少 800。这一步不要省，频点太稀的时候特征轨迹在拐弯处会被抹平，后面的净圈数统计直接出错。下一步要处理的排序问题，也是绘制环节里最容易翻车的地方。

3. 把开环传递矩阵离散成可绘制的特征轨迹：频点设计与排序算法

3.1 频率向量怎么取：log等距、覆盖范围与频点数

绘制广义奈奎斯特曲线时，频率向量必须按对数等距分布，不能用线性等距。系统穿越频率往往横跨几个十倍频程，线性分布会导致低频点数过密、高频点数稀疏，而高频段的相位特性又恰好决定轨迹怎么绕回原点。

经验取法是这样：先大致看对角通道的穿越频率范围，把频率下限设为穿越频率的 1/100 到 1/10，频率上限设为穿越频率的 10 倍以上。保守做法是直接取 1e-3 到 1e3 rad/s，点位数取 800 到 1200。每十倍频程至少有 100 个点，这个密度基本能保住轨迹的尖角。

# 推荐做法：对数等距频点，跨 6 个十倍频程，800 点 w_rad = np.geomspace(1e-3, 1e3, 800)

如果发现某段轨迹在关键拐弯处看起来像被“切了一刀”，不要只增加总点数，应该把频率区间收窄，只对穿越频率附近做局部加密。两个频段的重叠区域可以直接拼接，绘图时不影响曲线连续性。这个做法比无限加密省时间，也更容易判断频点是否够密。

注意：如果系统含有纯积分环节，最低频率不要取 0，取一个大于 0 的数值（比如 1e-3 甚至 1e-4），并在判据说明里单独记下虚轴极点，这部分在第 5 章会细讲。

3.2 模型来源两条路：解析传递函数与实测频率响应

有解析传递函数时，直接用第 2 章的方式逐频点代入即可。很多工程项目拿到的只有扫频数据，比如每个频率点上一组幅值和相角，没有拉普拉斯表达式。这时要先插值加密，再做逐点特征值分解。

实测频率响应插值有个常见的坑：用三次样条会在谐振峰附近过冲，幅频数据还好，相频数据一旦过冲，特征轨迹会出现不存在的局部圈，判据直接误报。我一般用保形插值（PCHIP）替代三次样条，它不会产生新的局部极值，在谐振峰附近更稳妥。

from scipy.interpolate import PchipInterpolator # f_exp: 实测扫频频率(Hz) # data_ij: 实测 G_ij(j*2*pi*f) 复数序列 def interp_G_entry(f_exp, data_ij, f_dense): re = PchipInterpolator(f_exp, data_ij.real)(f_dense) im = PchipInterpolator(f_exp, data_ij.imag)(f_dense) return re + 1j * im

注意插值对象是实部和虚部分别插值，不要直接对幅值和相角插值。相角数据在穿越 ±π 时会有跳变，不先做 unwrap 的话，插值出来的相角会在跳变沿产生畸变。如果扫频数据里相角已经呈现跳跃，先对原始相角做相位展开，再插值，最后转回复数。

插值完成后，每个频点都能得到完整的 2×2 复数矩阵，再进入特征值分解流程。实测数据的频点范围一般比解析模型窄，插值只能在你测量的频带内做，超出范围的特征轨迹没有意义，绘图时要明确标注有效频带。

3.3 特征值排序：不排序画出来就是一团乱麻

逐频点调用 numpy.linalg.eigvals 得到特征值后，频点 i 和频点 i+1 的特征值编号是独立的。理论上特征轨迹随 ω 连续变化，但数值计算返回的特征值顺序不保证相邻频点一致，直接把所有特征值依次连线，画出来就是多条轨迹相互交叉、断开的乱麻。

排序的思路是相邻频点做最优匹配：把频点 i 的特征值集合和频点 i+1 的特征值集合看成两组点，计算两两之间的距离，用二分图最小代价匹配求解。对小系统（2×2 或 3×3）速度非常快，对大系统也可以用同样的思路预筛后匹配。

from scipy.optimize import linear_sum_assignment def sort_trajectories(vals, w_rad): # vals: shape (n_freq, n_eig) # 返回按轨迹编号对齐的特征值数组 nf, nk = vals.shape order = np.zeros_like(vals, dtype=int) # 起点按实部降序，保证同一分支出发 order[0] = np.argsort(vals[0].real)[::-1] for i in range(1, nf): prev = vals[i-1, order[i-1]] # 上一频点已排序的特征值 cur = vals[i] # 当前频点原始特征值 cost = np.abs(prev[:, None] - cur[None, :])**2 # 距离平方成本矩阵 row, col = linear_sum_assignment(cost) # row=上帧编号 col=当前编号 perm = np.empty(nk, dtype=int) perm[col] = row # 当前 col 继承上帧 row 的轨迹号 order[i] = order[i-1][perm] return np.take_along_axis(vals, order, axis=1)

成本矩阵里 prev[:, None] 的形状是 (nk, 1)，cur[None, :] 是 (1, nk)，两者相减得到 nk×nk 的矩阵，每个元素表示上一帧某个特征值与当前帧某个特征值的距离平方。linear_sum_assignment 求解的是全局最小总代价，而不是每列独立找最近点，这样在两条轨迹距离接近的交叉区域不会出现来回跳号的情况。

排序算法输出后，还要做一步检查：在两条特征轨迹真正交叉的频点附近，特征值可能几乎重合，匹配存在歧义。这时候不是调排序算法，而是把频点加密，让交叉区域的两个特征值在数值上分开，排序才有稳定的依据。排序完成后，后面所有绘图和判据都基于排序后的数组，这比画完再去“修线”靠谱得多。

4. 把曲线画出来并完成稳定性判定：绘图代码与判读方法

4.1 完整绘制代码：补负频、标临界点、画方向

排序完成后，绘图前还有关键一步：补全负频支。前面说过实有理矩阵的特征值在负频率处取共轭，所以把正频支反转顺序再取共轭，拼接到正频支前面，就形成闭合曲线。只有闭合曲线才方便数净圈数，也才能直观看到轨迹怎么绕回起点。

import matplotlib.pyplot as plt def plot_characteristic_loci(vals_sorted, w_rad, kp): nf, nk = vals_sorted.shape trajectories = [] for k in range(nk): # 负频支：正频支反转顺序后取共轭 traj_neg = vals_sorted[::-1, k].conj() traj_pos = vals_sorted[:, k] traj_full = np.concatenate([traj_neg, traj_pos]) trajectories.append(traj_full) plt.plot(traj_full.real, traj_full.imag, lw=1.8, label=f"λ{k+1} (kp={kp})") # 起点画圆圈，终点画方块，区分方向 plt.plot(vals_sorted[0, k].real, vals_sorted[0, k].imag, 'o', ms=7) plt.plot(vals_sorted[-1, k].real, vals_sorted[-1, k].imag, 's', ms=7) plt.axhline(0, color='gray', lw=0.8) plt.axvline(0, color='gray', lw=0.8) plt.plot([-1], [0], 'rx', ms=10) # 临界点 plt.axis('equal') plt.grid(alpha=0.3) plt.legend() return trajectories

注意闭合顺序：traj_neg 是从最高频往最低频排列的共轭数据，它的第一个点在最高频处，最后一个点在 ω=0 处；traj_pos 的第一个点也在 ω=0 处。两个数组在 ω=0 处的点重合，整个轨迹从 ω=−∞ 扫到 ω=+∞。matplotlib 里 axis('equal') 很重要，坐标轴不等比时，轨迹绕临界点的视觉距离会被拉伸，容易误判裕量。

起点和终点用不同 marker，是为了在轨迹密集时快速分清扫频方向。如果还想加方向箭头，可以在轨迹每间隔 100 个点左右画一个小箭头，但箭头太密反而遮挡轨迹，我一般只保留起点点和终点方块，配合图例里的 ω 递增方向说明就够了。

4.2 如何数净环绕次数：补镜像后比圈数

闭合曲线数圈数，最稳妥的方式是对轨迹相对 (−1, j0) 的角度做 unwrap，看首尾角度变化量除以 2π 是多少。因为闭合格局下，角度从终点回到起点时，差值一定是 2π 的整数倍。

def encirclement_net(traj): # traj: 补全后的闭合特征轨迹，逆时针环绕为正 z = traj + 1 # 相对临界点的复数坐标 ang = np.unwrap(np.angle(z)) return int(np.round((ang[-1] - ang[0]) / (2 * np.pi)))

np.unwrap 会自动累加角度跳变。比如轨迹绕着 (−1, j0) 逆时针走了一圈，角度从 0 附近逐渐增大到 2π，unwrap 后终值约等于初值加 6.283，除以 2π 就是 1。顺时针则得到 −1。注意如果频点太稀，轨迹在尖角处相位变化超过 π，unwrap 会发生误判，这就是为什么前面反复强调频点数不能省。

判据落地很简单：把所有特征轨迹的净圈数加起来得到 N，再看开环右半平面极点 P，Z = P + N 为 0 则闭环稳定。以第 2 章那个 2×2 对象为例，取 kp=0.3 时两条特征轨迹都远离 (−1, j0)，N = 0，P = 0，闭环稳定；把 kp 逐步增大到 0.8 附近，轨迹开始逼近临界区域，如果现场能观察到发散，多半是某条轨迹已经完成了绕圈。

4.3 判读流程与绘图参数

判读不要只盯着图看，按固定顺序走，能少很多误判。

首先计算开环右半平面极点数 P，用解析方式求 L(s) 的全部极点，挑实部大于 0 的；不要想当然认为被控对象稳定 P 就是 0。其次逐条统计特征轨迹的净圈数，加总得到 N。最后算 Z = P + N，Z=0 说明闭环稳定，否则右半平面极点数就是 Z 个。判读完成后，再做一次时域阶跃仿真验证，两条路径的结果对得上才算数。

曲线多的时候，建议每画完一版就保存一个 PNG 快照，标注 kp 值。后面调参数做对比时，快照比代码里反复改值更直观。这是我做多变量频域分析的习惯，也是排查问题时的第一手资料。

5. 广义奈奎斯特曲线绘制避坑：五个典型翻车现场与排查顺序

5.1 频点太稀：特征轨迹的“尖角”被抹平，包围次数数错

现象：画出来的特征轨迹在高频段或者穿越频带附近显得过于平滑，和对象实际频率特性对不上；净圈数统计结果在 kp 变化时跳变不连续。

原因：频率向量点位数不足，或者用了线性等距频率分布。特征轨迹在特征值接近趋近的频段会有快速转向，采样点跨过这个区域时，转向圆弧被一条直线替代，unwrap 统计圈数时直接丢了一圈。

解决：把频率向量改成对数等距，并保证每十倍频程至少 100 个点。如果局部还不够，就在转向密集区间单独加密，比如 0.1 到 10 rad/s 之间再插入一组更密的频点，和全局频点合并后重画。

5.2 特征值编号跳变：曲线断裂、交叉，负频补全出现镜像错位

现象：同一颜色的轨迹在某处突然跳到另一条轨迹上，曲线呈锯齿状；补负频镜像后，闭合曲线围绕的区域明显不对称。

原因：每个频点独立调用 eigvals，特征值的返回顺序不保证连续，未做排序就直接连线。特征值在交叉频点附近数值相近，数值库返回顺序随机，轨迹编号自然就跳了。

解决：先把第 3 章的排序算法接到特征值分解后面，对每个频点做相邻匹配。排序后最直观的检查方式是看图：每条轨迹颜色连续，交叉处两条不同颜色轨迹可以靠近但不会互换颜色。

5.3 轨迹正好穿过临界点(−1, j0)：判据退化

现象：轨迹精确经过 (−1, j0)，角度 unwrap 产生奇异，净圈数统计结果不确定；kp 微调一点，判定结果就从稳定翻到不稳定。

原因：特征轨迹与临界点重合意味着 det(I + L(jω)) = 0，闭环极点在虚轴上，系统处于临界稳定状态。这不是绘图错误，而是系统本身就处在边界。频点离散化可能掩盖真正的接触点，导致看起来像是“擦边”。

解决：先判定是否真的穿过，加密局部频率，把接触点找出来；然后微调 kp（比如 ±2%），看轨迹向临界点左侧还是右侧移动。如果轨迹正好在 kp=0.73 穿过临界点，把 kp 调成 0.72 或 0.74 再画一版，判据的边际方向就清楚了。

5.4 开环含积分器：虚轴极点让判据失效

现象：频率接近 0 时，特征轨迹向实轴无穷远处冲出去，曲线无法闭合，unwrap 结果完全没有规律。

原因：积分环节把 L(s) 的极点放在 s=0，正好在 Nyquist 围线上。经典奈奎斯特判据要求围线绕开虚轴极点，用半径为 ε 的小半圆绕行，映射到曲线里会出现一个半径无穷大的圆弧。逐点绘制到 ω=0 时，这个操作被忽略了。

解决：工程上最简单的做法是频率下限取一个大于 0 的数值，比如 1e-3 rad/s，并单独记录虚轴极点个数 q。判据需要修正：稳定条件从 Z = P + N 变成 Z = P + N + q/2 的对应修正项，不同书上的约定不一样，用之前先明确自己的符号体系。我习惯的做法是把积分环节先提出来，单独分析完再合并结果，避免在判据里混用。

5.5 开环右半平面极点P读错：判据和仿真对不上

现象：广义奈奎斯特判据判定闭环稳定，但时域 step 仿真发散；或者判据判定不稳定，现场却稳定运行。

原因：P 的统计漏项。最常见的是只统计了被控对象 G(s) 的右半平面极点，忘了控制器 K(s) 里也可能有；另一个常见问题是把虚轴上的极点也算进 P，符号体系混乱。

解决：直接用数值方法求 L(s) = G(s)K(s) 的传递函数极点并筛选实部大于 0 的项，不要手数。把 P 的统计过程和特征轨迹图画在同一张图里，作为输出数据一份存档，调参时每次改 K(s) 都要重算 P，不要复用旧值。记住 GNC 里 P 是开环回路整体的极点，不是被控对象单独的极点。

6. 进阶：用特征轨迹定位临界频率，再做时域闭环验证

6.1 从穿越点反解临界频率

广义奈奎斯特曲线不仅能判定稳定，还能定位控制器最需要补偿的频率点。轨迹与负实轴交点（实部小于 −1）对应的频率就是潜在的临界频率，是相位裕量和增益裕量的直接体现。

def locate_crossing(w_rad, traj): # 查找特征轨迹与负实轴左侧的交点频率 for k in range(traj.shape[1]): y = traj[:, k].imag x = traj[:, k].real for i in range(len(w_rad) - 1): if y[i] * y[i+1] < 0 and x[i] < -1 and x[i+1] < -1: w_c = w_rad[i] + (w_rad[i+1] - w_rad[i]) * (-y[i]) / (y[i+1] - y[i]) print(f"λ{k+1} 在 {w_c:.3f} rad/s 穿过负实轴左侧")

这个频率信息比单纯看曲线更有用。比如定位到临界频率在 0.8 rad/s，而耦合通道的转折频率恰好也在 0.8 rad/s，那就说明耦合支路主导了稳定性边界，控制器设计应该优先降耦合增益或者在该频段加超前补偿，而不是盲目调主对角增益。

6.2 把判据结论送进时域仿真验证

频域判据说到底是个数学结论，工程验收还得看时域响应。我每次都把 GNC 判据和闭环阶跃仿真并排跑一遍，两个结论一致才敢继续调参数。

# 伪代码：构造闭环状态空间后做阶跃仿真验证 # from scipy.signal import lsim # sys_cl = ... 由 G(s) 和 K(s) 组成的闭环状态空间 # t = np.linspace(0, 20, 2000) # y, t, x = lsim(sys_cl, U=step_input, T=t) # 观察两个输出通道是否收敛，是否出现等幅振荡

如果判据说稳定但仿真发散，优先检查 P 是否统计错了，再看是否漏掉了虚轴极点的修正项；如果判据说稳定且仿真也稳定，再看裕量够不够。特征轨迹离临界点越远，时域响应阻尼越好，这个距离就是多回路控制器最直观的稳定性余量。

这套方法我用了很多年，最大的教训是别跳过排序和负频补全直接看图像，也别在频点数上省时间——多数“画出来不对劲”的情况，最后都能归到这两件事上。先把闭环曲线画闭合、把编号对齐，再去读稳定性，后面的判断就会顺很多。希望帮到你。

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