Python实战：用control库5步搞定根轨迹绘制（附避坑指南）

如果你正在学习自动控制原理，或者从事控制系统的设计与调试工作，那么“根轨迹”这个概念对你来说一定不陌生。它就像一张动态地图，清晰地展示了当系统某个参数（通常是开环增益）变化时，闭环极点在复平面上移动的轨迹。传统教材里，我们往往需要背诵十条绘制法则，在纸上进行繁琐的几何作图。但现在，我们有了更强大的工具——Python。

这篇文章就是为你准备的，无论你是自动化专业的学生，还是需要快速验证设计方案的工程师。我们将彻底抛开枯燥的理论推导，聚焦于如何用Python的control库，在短短五步之内，从零生成一张专业、精确的根轨迹图。更重要的是，我会结合自己踩过的坑，分享那些教科书上不会写的实战技巧和避坑指南，比如如何处理复数运算、为什么你的渐近线画不出来、以及如何从根轨迹中直接“读”出系统的动态性能。让我们告别手工作图，用代码来驾驭控制系统的分析与设计。

1. 环境搭建与库安装：迈出坚实的第一步

在开始绘制根轨迹之前，我们需要一个合适的“工作台”。对于控制系统的分析与仿真，Python生态中的control库是不二之选。它提供了类似MATLAB Control System Toolbox的接口，但完全免费且开源。不过，它的安装和配置有一些小细节需要注意，处理不好可能会导致后续步骤无法进行。

首先，确保你有一个Python环境（推荐3.8及以上版本）。打开你的终端或命令提示符，最直接的安装命令是使用pip：

pip install control

但是，这里有一个常见的坑：control库的核心功能依赖于numpy和scipy进行数值计算，而matplotlib用于绘图。如果你的Python环境是全新的，或者之前安装的scipy版本不兼容，可能会报错。一个更稳妥的做法是使用conda（如果你使用Anaconda发行版）或者创建一个干净的虚拟环境。

# 使用conda安装（推荐用于科学计算环境）
conda install -c conda-forge control

# 或者，使用pip并指定依赖版本以确保兼容性
pip install numpy scipy matplotlib
pip install control

安装完成后，不要急着写代码。先进行一个简单的“冒烟测试”，验证库是否正常工作，并解决一个中文用户经常遇到的问题——图形中的字体显示。

import control
import matplotlib.pyplot as plt
import numpy as np

# 测试1：检查control库版本和基本功能
print(f"control库版本: {control.__version__}")
# 尝试创建一个简单的传递函数
try:
    sys_test = control.TransferFunction([1], [1, 2, 1])
    print("传递函数创建成功:", sys_test)
except Exception as e:
    print("创建传递函数失败:", e)

# 测试2：配置matplotlib以正确显示中文和负号
# 这是避免图形出现乱码或框框的关键步骤
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']  # 设定中文字体，备选英文字体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示为方块的问题
print(" matplotlib 字体配置完成。")

如果以上代码能顺利运行并打印出版本信息，那么你的基础环境就准备好了。我建议将字体配置代码保存为一个单独的模块或放在脚本开头，这样所有绘图都能保持一致。

2. 定义系统模型：从传递函数到Python对象

根轨迹描绘的是闭环极点随开环增益变化的轨迹，因此一切分析的起点都是系统的开环传递函数。在control库中，我们使用control.TransferFunction对象来表示它。这一步看似简单，但输入格式的细微错误会导致完全错误的结果。

传递函数通常表示为 \( G(s) = \frac{N(s)}{D(s)} \)，其中 \( N(s) \) 和 \( D(s) \) 是s的多项式。在Python中，我们用系数列表来定义这些多项式，列表的顺序是从最高次幂到常数项。这是最容易出错的地方。

假设我们要分析一个经典的三阶系统：\( G(s) = \frac{K}{s(s+1)(s+5)} \)。首先需要将其展开为多项式形式： \( G(s) = \frac{K}{s^3 + 6s^2 + 5s} \)。 因此，分子多项式是 \( K \cdot 1 \)，分母多项式是 \( s^3 + 6s^2 + 5s + 0 \)。

对应的Python代码如下：

# 定义分子和分母多项式的系数
# 分子: 1 (代表K*1，K在绘制根轨迹时会自动变化)
numerator = [1]  # 等价于 1*s^0
# 分母: s^3 + 6s^2 + 5s + 0
denominator = [1, 6, 5, 0]  # 系数依次为: s^3, s^2, s^1, s^0

# 创建传递函数对象
sys = control.TransferFunction(numerator, denominator)
print("开环传递函数 G(s) = ", sys)

运行后，你应该能看到类似 1/(s^3 + 6 s^2 + 5 s) 的输出。这就表示对象创建成功了。

为了应对更复杂的系统，这里有一个实用技巧表格，总结了不同形式传递函数的定义方法：

传递函数形式	数学表达式	分子系数 (num)	分母系数 (den)	说明
增益+极点	\( \frac{K}{s+a} \)	[K]	[1, a]	注意分母是s+a
二阶系统	\( \frac{K}{s^2 + 2\zeta\omega_n s + \omega_n^2} \)	[K]	[1, 2*zeta*wn, wn**2]	常用于震荡系统分析
含零点系统	\( \frac{K(s+z)}{s(s+p)} \)	[K, K*z]	[1, p, 0]	分子也要展开
零极点形式	\( K\frac{(s+2)}{(s+1)(s+3)} \)	[K, 2*K]	[1, 4, 3]	先乘开(s+1)(s+3)=s^2+4s+3

注意：在定义分母时，尤其要留意常数项。例如，对于 \( s(s+1) \)，展开是 \( s^2 + s \)，常数项为0，因此分母是 [1, 1, 0]。如果错误地写成 [1, 1]，就相当于丢掉了s这个因子，系统阶数会出错，导致根轨迹分支数减少。

创建好系统对象后，我们可以先用control.pole()和control.zero()函数查看其开环零极点，这对后续理解根轨迹的起点和终点很有帮助。

poles = control.pole(sys)
zeros = control.zero(sys)
print(f"开环极点: {poles}")
print(f"开环零点: {zeros}")

对于我们的例子，应该会看到极点在 [0, -5, -1]，零点列表为空。这符合预期：根轨迹将从这三个点出发。

3. 核心绘制与基本解读：生成你的第一张根轨迹图

有了系统对象，绘制根轨迹就变得异常简单。control.rlocus()函数是完成这项工作的核心。但是，直接调用control.rlocus(sys)并坐等出图，可能会让你错过很多细节，并且在需要进一步分析时束手无策。

让我们先看看最基础的绘制方法：

# 方法1：最简单的方式，自动生成绘图
plt.figure(figsize=(10, 8))
control.rlocus(sys)
plt.title('系统根轨迹图 - 基础绘制')
plt.grid(True)
plt.show()

执行这段代码，你会得到一张标准的根轨迹图。x轴代表实部，y轴代表虚部。图上应该有三条曲线（对应三个极点），它们从开环极点（0, -1, -5）开始，随着K增大而移动。其中两条轨迹最终穿过虚轴进入右半平面，这预示着系统在K值超过某个临界值后会变得不稳定。

然而，这种“黑箱”式的绘图往往不能满足我们的需求。比如，我们想知道特定K值对应的闭环极点是什么，或者想自定义图形的样式。这时，就需要使用control.rlocus的“非绘图”模式，并获取其返回值。

# 方法2：高级方式，获取数据并自定义绘图
plt.figure(figsize=(10, 8))
# 设置计算增益K的范围，避免无穷大导致的计算问题
kvect = np.logspace(-2, 3, 1000)  # 从10^-2到10^3，取1000个对数间隔点
# rlocus返回两个值：根轨迹点数组 和 对应的增益数组
roots, gains = control.rlocus(sys, kvect=kvect, plot=False)

# 现在我们可以完全控制如何绘图
for i in range(roots.shape[1]):  # 遍历每一支根轨迹
    plt.plot(roots[:, i].real, roots[:, i].imag, 'b-', linewidth=1.5, label='根轨迹' if i==0 else "")
    # 标记起点（K=0时的极点）
    plt.plot(roots[0, i].real, roots[0, i].imag, 'ko', markersize=8)
    # 标记终点（K很大时的趋向）
    plt.plot(roots[-1, i].real, roots[-1, i].imag, 'k^', markersize=8)

# 添加坐标轴和网格
plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)  # 实轴
plt.axvline(x=0, color='k', linestyle='-', alpha=0.3)  # 虚轴
plt.xlabel('实部 (σ)')
plt.ylabel('虚部 (jω)')
plt.title('自定义系统根轨迹图（含起止点标记）')
plt.grid(True, alpha=0.3)
plt.axis('equal')  # 保证x轴和y轴比例相同，圆形才不会变椭圆
plt.legend()
plt.show()

这段代码给了我们更大的灵活性。roots变量是一个三维数组，其形状大致为 (len(kvect), n)，其中n是系统阶数。roots[:, i] 就是第i支根轨迹上所有点的复数坐标。

现在，我们可以轻松回答“当K=10时，闭环极点在哪儿？”这个问题：

# 寻找最接近K=10的增益索引
target_K = 10
idx = np.argmin(np.abs(gains - target_K))
poles_at_K10 = roots[idx, :]
print(f"当 K ≈ {gains[idx]:.2f} 时，闭环极点为：")
for p in poles_at_K10:
    print(f"  {p:.4f}")

通过这种方式，根轨迹从一张静态的图片，变成了一个可以交互查询的数据集。你可以修改kvect的范围和密度，来重点关注你感兴趣的增益区间（例如，临界稳定点附近）。

4. 高级分析与常见陷阱：超越基础绘图

仅仅画出轨迹还不够，控制工程师需要从图中提取关键信息：分离点、与虚轴的交点（临界稳定增益）、渐近线方向等。这些是评估系统稳定裕度和动态性能的核心。手动计算这些点非常繁琐，但用Python辅助，我们可以半自动化地完成。

4.1 计算分离点/汇合点

分离点是根轨迹在实轴上离开或汇合的点，对应闭环特征方程的重根。我们可以通过求解方程 \( dK/ds = 0 \) 来找到它们。对于系统 \( 1 + KG(s) = 0 \)，有 \( K = -1/G(s) \)。分离点满足 \( dK/ds = d(-1/G(s))/ds = 0 \)。

import sympy as sp

# 使用符号计算寻找分离点
s = sp.symbols('s')
# 定义开环传递函数 G(s) = 1/(s*(s+1)*(s+5))
G_s = 1 / (s * (s + 1) * (s + 5))
# K = -1/G(s)
K_expr = -1 / G_s
# 求导并解方程
dK_ds = sp.diff(K_expr, s)
solutions = sp.solve(dK_ds, s)
print("可能的分离点（符号解）: ", solutions)

# 过滤出实数解，并且该点在实轴的根轨迹段上
real_solutions = [complex(sol.evalf()) for sol in solutions if sp.im(sol) == 0]
print("实数解: ", real_solutions)

# 验证：这些点是否在根轨迹上（是否满足相角条件）
# 一个简易验证：检查该点右侧实轴上的零极点个数是否为奇数
for sol in real_solutions:
    sigma = sol.real
    # 对于本例，实轴上的根轨迹区间是(-∞, -5] 和 [-1, 0]
    if (-5 >= sigma) or (-1 <= sigma <= 0):
        print(f"σ = {sigma:.3f} 是有效的分离点。")
        # 计算对应的K值
        K_val = abs(1/float(G_s.subs(s, sigma)))
        print(f"  对应的增益 K ≈ {K_val:.3f}")

4.2 计算与虚轴的交点（临界稳定点）

根轨迹与虚轴的交点意味着系统处于临界稳定状态，这时的增益 \( K_{crit} \) 是稳定性的边界。我们可以通过令闭环特征方程的 \( s = j\omega \) 来求解。

# 系统特征方程: 1 + K*G(s) = 0 -> s(s+1)(s+5) + K = 0
# 即 s^3 + 6s^2 + 5s + K = 0
# 令 s = jω，代入： (jω)^3 + 6(jω)^2 + 5(jω) + K = 0
# 整理： (-jω^3) - 6ω^2 + 5jω + K = 0
# 实部： -6ω^2 + K = 0
# 虚部： -ω^3 + 5ω = 0 -> ω(5 - ω^2) = 0

# 解方程
omega_solutions = []
# 从虚部方程解ω
# ω = 0 是一个解（对应起点），通常我们关心非零解
omega_nonzero = np.sqrt(5)  # ω^2 = 5
omega_solutions.append(omega_nonzero)
omega_solutions.append(-omega_nonzero)  # 共轭，通常只取正数分析

for omega in omega_solutions:
    if omega > 0:
        # 代入实部方程求K
        K_crit = 6 * omega**2
        print(f"根轨迹与虚轴交点: s = ±j{omega:.3f}")
        print(f"  对应的临界增益 K_crit = {K_crit:.3f}")
        # 在图中标记这个点
        plt.plot(0, omega, 'r*', markersize=15, label='临界稳定点')
        plt.plot(0, -omega, 'r*', markersize=15)

4.3 绘制渐近线

当系统阶数n大于零点数m时，有n-m条根轨迹趋向于无穷远，其渐近线提供了轨迹走向的大致趋势。渐近线与实轴的交点 \( \sigma_a \) 和夹角 \( \theta_a \) 有固定公式。很多初学者发现control库不自动绘制渐近线，我们可以手动补上。

# 计算渐近线参数
n = len(poles)  # 极点数
m = len(zeros)  # 零点数
num_asymp = n - m

if num_asymp > 0:
    # 交点 sigma_a = (极点之和 - 零点之和) / (n - m)
    sigma_a = (np.sum(poles) - np.sum(zeros)) / num_asymp
    print(f"渐近线与实轴交点: σ_a = {sigma_a:.3f}")
    
    # 夹角 theta_a = (2k+1)*180° / (n-m), k=0,1,..., (n-m-1)
    angles_deg = []
    for k in range(num_asymp):
        angle = (2*k + 1) * 180.0 / num_asymp
        # 通常将角度规范到 (-180°, 180°] 区间
        if angle > 180:
            angle -= 360
        angles_deg.append(angle)
    print(f"渐近线夹角: {angles_deg}")
    
    # 在图上绘制渐近线（用虚线表示）
    plt.axvline(x=sigma_a, color='gray', linestyle=':', alpha=0.5, label='渐近线交点')
    # 绘制每条渐近线（取一段线段示意）
    r = 10  # 线段长度
    for angle in angles_deg:
        angle_rad = np.deg2rad(angle)
        dx = r * np.cos(angle_rad)
        dy = r * np.sin(angle_rad)
        plt.plot([sigma_a, sigma_a + dx], [0, dy], 'g--', alpha=0.5, linewidth=1)
        plt.plot([sigma_a, sigma_a + dx], [0, -dy], 'g--', alpha=0.5, linewidth=1)

将这些分析元素整合到一张图中，你就能得到一张信息量远超默认输出的专业级根轨迹分析图。

5. 从根轨迹到系统性能：建立直观联系

绘制根轨迹的最终目的，是为了分析和设计系统。我们需要建立“图上位置”与“时域性能”之间的直观联系。这主要通过主导极点的概念来实现。主导极点是指离虚轴最近、且附近没有闭环零点的极点，它们主导了系统的瞬态响应。

对于二阶系统，极点位置 \( s = -\zeta\omega_n \pm j\omega_n\sqrt{1-\zeta^2} \) 与阻尼比 \( \zeta \)、自然频率 \( \omega_n \)、超调量、调节时间有精确的关系。对于高阶系统，如果存在一对共轭主导极点，我们可以近似用二阶系统的公式来估算性能。

让我们写一段代码，在根轨迹图上画出等阻尼比线（射线）和等自然频率线（圆弧），并演示如何根据给定的极点位置估算性能指标。

# 在根轨迹图上添加等阻尼比线和等自然频率线
zeta_values = [0.2, 0.5, 0.707, 0.9]  # 常见的阻尼比值
omega_n = 3  # 选择一个自然频率值用于画圆弧

plt.figure(figsize=(12, 10))
# 先绘制根轨迹
roots, gains = control.rlocus(sys, kvect=np.logspace(-2, 3, 500), plot=False)
for i in range(roots.shape[1]):
    plt.plot(roots[:, i].real, roots[:, i].imag, 'b-', alpha=0.7)

# 1. 绘制等阻尼比线 (射线)
for zeta in zeta_values:
    theta = np.arccos(zeta)  # 阻尼比对应的角度
    # 画一条从原点出发的射线
    x_line = np.linspace(-6, 0, 100)
    y_line_pos = x_line * np.tan(theta)  # 上半平面
    y_line_neg = -y_line_pos  # 下半平面
    # 只画左半平面部分
    valid_idx = x_line < 0
    plt.plot(x_line[valid_idx], y_line_pos[valid_idx], 'r--', alpha=0.5, linewidth=0.8)
    plt.plot(x_line[valid_idx], y_line_neg[valid_idx], 'r--', alpha=0.5, linewidth=0.8)
    # 添加阻尼比标签
    label_x = -2
    label_y = label_x * np.tan(theta) + 0.2
    plt.text(label_x, label_y, f'ζ={zeta}', fontsize=9, color='red', alpha=0.8)

# 2. 绘制等自然频率线 (圆弧， s^2 + 2ζω_n s + ω_n^2 = 0 的根轨迹是圆)
theta_circle = np.linspace(0, np.pi, 200)
# 圆的参数方程：实部 = -ζω_n, 虚部 = ω_n sinθ， 但这里我们画一个以原点为圆心，半径为ω_n的圆弧在左半平面
# 更准确的是画等ω_n线，即极点与原点的距离为ω_n
circle_x = -omega_n * np.cos(theta_circle)  # 确保在左半平面
circle_y = omega_n * np.sin(theta_circle)
plt.plot(circle_x, circle_y, 'm--', alpha=0.5, linewidth=0.8, label=f'ω_n ≈ {omega_n}')
plt.text(-omega_n*0.8, omega_n*0.7, f'ω_n={omega_n}', fontsize=9, color='purple', alpha=0.8)

plt.axhline(0, color='k', linestyle='-', alpha=0.2)
plt.axvline(0, color='k', linestyle='-', alpha=0.2)
plt.xlabel('实部 (σ)')
plt.ylabel('虚部 (jω)')
plt.title('根轨迹图与等阻尼比、等自然频率线')
plt.grid(True, alpha=0.3)
plt.axis('equal')
plt.legend(['根轨迹', '等阻尼比线', '等自然频率线'])
plt.show()

现在，假设我们从根轨迹上选择了一个点，例如通过交互式操作或计算得到一个闭环极点 s = -0.5 + 1.5j。我们可以快速估算其对应的性能：

# 性能估算函数
def estimate_performance_from_pole(pole):
    """从单个复数极点（假设是主导共轭极点中的一个）估算二阶系统性能"""
    sigma = abs(pole.real)
    omega_d = abs(pole.imag)
    omega_n = abs(pole)  # 模长
    zeta = sigma / omega_n
    
    # 计算超调量 (Percent Overshoot, PO)
    if 0 < zeta < 1:
        PO = np.exp(-zeta * np.pi / np.sqrt(1 - zeta**2)) * 100
    else:
        PO = 0 if zeta >= 1 else np.nan
    
    # 计算调节时间 (2% criterion, Ts)
    Ts = 4.6 / sigma  # 当 zeta < 0.7 时比较准确
    
    # 计算峰值时间 (Tp)
    if omega_d > 0:
        Tp = np.pi / omega_d
    else:
        Tp = np.nan
        
    return zeta, omega_n, PO, Ts, Tp

# 示例：估算极点 s = -0.5 + 1.5j 的性能
example_pole = complex(-0.5, 1.5)
zeta_est, wn_est, PO_est, Ts_est, Tp_est = estimate_performance_from_pole(example_pole)

print("=== 从极点位置估算系统性能 ===")
print(f"极点: {example_pole:.3f}")
print(f"估算阻尼比 ζ: {zeta_est:.3f}")
print(f"估算自然频率 ω_n: {wn_est:.3f} rad/s")
print(f"估算超调量 PO: {PO_est:.1f}%")
print(f"估算调节时间 (2%) Ts: {Ts_est:.3f} s")
print(f"估算峰值时间 Tp: {Tp_est:.3f} s")

通过这样的分析，根轨迹就从抽象的数学曲线，变成了一个强大的设计工具。你可以根据期望的超调量（比如小于20%）和调节时间，在图上反推出主导极点应该位于哪个区域（例如，位于ζ=0.5的等阻尼比线附近），进而确定所需的开环增益K的范围。

最后，别忘了结合时域仿真来验证你的分析。control库可以方便地计算阶跃响应：

# 验证：选择一个增益K，查看闭环系统的阶跃响应
K_design = 10  # 假设我们选择K=10
closed_loop_sys = control.feedback(K_design * sys, 1)  # 单位负反馈闭环
t, y = control.step_response(closed_loop_sys, T=np.linspace(0, 20, 1000))

plt.figure(figsize=(10, 5))
plt.plot(t, y)
plt.title(f'闭环系统阶跃响应 (K={K_design})')
plt.xlabel('时间 [s]')
plt.ylabel('输出')
plt.grid(True)
plt.show()

将根轨迹分析、性能估算和时域仿真三者结合，你就能对系统的行为有一个全面而深刻的理解。在实际项目中，我通常会快速绘制根轨迹，根据性能要求确定一个大致的增益范围，然后在这个范围内精细调整，并用时域响应做最终验证。这种工作流极大地提高了设计效率，也减少了对直觉和经验的过度依赖。