别再死记硬背了！用Python+Matplotlib手把手教你动态绘制根轨迹（附完整代码）-程序员充电站

用Python+Matplotlib动态绘制根轨迹：从零实现自控原理可视化工具

打开你的Python环境，我们直接进入正题。传统自控原理教材中那些晦涩的根轨迹绘制法则，今天将通过代码的逐行实现变得触手可及。这不是简单的理论复述，而是一个完整的交互式可视化项目——你将亲手构建一个能动态响应参数变化的根轨迹绘制器，在代码运行过程中直观理解极点移动、系统稳定性等核心概念。

1. 环境配置与基础理论准备

在开始编码前，确保你的Python环境已安装以下库：

pip install numpy matplotlib scipy control

根轨迹的本质是描述系统闭环极点随开环增益K变化的运动轨迹。考虑一个典型二阶系统：

$$ G(s) = \frac{K(s+z)}{(s+p_1)(s+p_2)} $$

当K从0→∞变化时，闭环特征方程1+G(s)H(s)=0的根在复平面上的移动轨迹就是我们要绘制的根轨迹。这种可视化呈现比纯数学推导更能直观展示：

稳定性边界：轨迹是否跨越虚轴
动态响应特性：极点位置与阻尼比、自然频率的关系
参数敏感性：零点/极点位置微小变化对系统的影响

提示：本文所有代码块都可直接复制到Jupyter Notebook中分段执行，建议边阅读边实践

2. 核心算法实现：从传递函数到轨迹计算

我们首先构建计算根轨迹的底层引擎。创建一个RootLocusPlotter类来封装所有核心功能：

import numpy as np from scipy import signal import matplotlib.pyplot as plt from matplotlib.widgets import Slider class RootLocusPlotter: def __init__(self, num, den): """初始化传递函数系数 num: 分子多项式系数列表，如[1,2]表示s+2 den: 分母多项式系数列表，如[1,3,2]表示s²+3s+2 """ self.num = np.array(num) self.den = np.array(den) self.k_values = np.logspace(-2, 4, 500) # K值范围：0.01到10000 self.fig, self.ax = plt.subplots(figsize=(10, 8)) def compute_root_locus(self): """计算根轨迹数据""" roots = [] for k in self.k_values: closed_loop_den = self.den + k * np.pad(self.num, (len(self.den)-len(self.num), 0), 'constant') root = np.roots(closed_loop_den) roots.append(root) return np.array(roots)

这个基础版本已经能处理常规根轨迹计算。测试一个简单系统：

# 示例：G(s)=K(s+1)/(s²+3s+2) rl_plotter = RootLocusPlotter(num=[1, 1], den=[1, 3, 2]) roots = rl_plotter.compute_root_locus() # 可视化结果 for i in range(roots.shape[1]): plt.plot(roots[:, i].real, roots[:, i].imag, 'b-') plt.xlabel('Real Axis') plt.ylabel('Imaginary Axis') plt.grid(True) plt.show()

3. 动态可视化：让根轨迹"活"起来

静态图像无法展现参数变化的动态过程，我们引入Matplotlib的交互功能：

def plot_interactive_rl(num, den): """创建交互式根轨迹图""" fig, ax = plt.subplots(figsize=(10, 8)) plt.subplots_adjust(bottom=0.25) # 初始绘制 rl = RootLocusPlotter(num, den) roots = rl.compute_root_locus() lines = [] for i in range(roots.shape[1]): line, = ax.plot(roots[:, i].real, roots[:, i].imag, 'b-') lines.append(line) # 添加零点/极点标记 zeros = np.roots(num) poles = np.roots(den) ax.scatter(zeros.real, zeros.imag, marker='o', facecolors='none', edgecolors='r', s=100) ax.scatter(poles.real, poles.imag, marker='x', color='r', s=100) # 创建K值滑块 ax_k = plt.axes([0.25, 0.1, 0.65, 0.03]) k_slider = Slider(ax_k, 'K值', 0.01, 100, valinit=1, valfmt='%0.2f') def update(val): k = k_slider.val closed_loop_den = den + k * np.pad(num, (len(den)-len(num), 0), 'constant') current_roots = np.roots(closed_loop_den) # 更新当前极点位置标记 if hasattr(update, 'current_markers'): for m in update.current_markers: m.remove() update.current_markers = [ ax.scatter(r.real, r.imag, color='g', s=80) for r in current_roots] fig.canvas.draw_idle() k_slider.on_changed(update) ax.set_xlabel('Real Axis') ax.set_ylabel('Imaginary Axis') ax.grid(True) plt.show() # 使用示例 plot_interactive_rl(num=[1, 2], den=[1, 5, 6, 0])

现在拖动滑块，你将看到：

绿色圆点实时显示当前K值对应的闭环极点位置
红色×和○分别表示开环极点和零点
蓝色曲线展示完整的根轨迹

4. 高级特性实现：分离点与渐近线计算

专业级的根轨迹工具需要包含关键特征点的计算。我们扩展RootLocusPlotter类：

class AdvancedRootLocusPlotter(RootLocusPlotter): def find_breakaway_points(self): """计算分离/汇合点""" # 构造特征多项式的导数方程 A = np.polymul(self.num, np.polyder(self.den)) B = np.polymul(self.den, np.polyder(self.num)) break_eq = np.polysub(A, B) possible_points = np.roots(break_eq) # 筛选实数且在根轨迹上的点 valid_points = [] for p in possible_points: if np.isreal(p): p = p.real # 验证是否在根轨迹上（实轴右侧零极点数为奇数） right_zeros = sum(z.real >= p for z in np.roots(self.num) if np.isreal(z)) right_poles = sum(pole.real >= p for pole in np.roots(self.den) if np.isreal(pole)) if (right_zeros - right_poles) % 2 == 1: valid_points.append(p) return valid_points def compute_asymptotes(self): """计算渐近线参数""" m = len(np.roots(self.num)) n = len(np.roots(self.den)) if n <= m: return None # 渐近线交点 sigma = (sum(np.roots(self.den)) - sum(np.roots(self.num))) / (n - m) # 渐近线角度 angles = [(2*k+1)*np.pi/(n-m) for k in range(n-m)] return sigma, angles

将这些特征添加到可视化中：

def plot_with_features(num, den): rl = AdvancedRootLocusPlotter(num, den) roots = rl.compute_root_locus() fig, ax = plt.subplots(figsize=(10, 8)) # 绘制基本轨迹 for i in range(roots.shape[1]): ax.plot(roots[:, i].real, roots[:, i].imag, 'b-') # 标记分离点 break_points = rl.find_breakaway_points() for bp in break_points: ax.scatter(bp, 0, color='purple', s=100, marker='s') # 绘制渐近线 asymptotes = rl.compute_asymptotes() if asymptotes: sigma, angles = asymptotes for angle in angles: x = np.linspace(sigma, sigma+10*np.cos(angle), 100) y = 10*np.sin(angle) * (x - sigma) / 10 ax.plot(x, y, 'r--', alpha=0.5) # 标记零极点 zeros = np.roots(num) poles = np.roots(den) ax.scatter(zeros.real, zeros.imag, marker='o', facecolors='none', edgecolors='r', s=100) ax.scatter(poles.real, poles.imag, marker='x', color='r', s=100) ax.set_xlabel('Real Axis') ax.set_ylabel('Imaginary Axis') ax.grid(True) plt.show() # 示例：展示所有高级特征 plot_with_features(num=[1], den=[1, 3, 2, 0])

5. 实战应用：系统性能分析与设计

通过根轨迹可以直接评估系统性能。我们创建一个分析工具：

def system_performance_analysis(num, den): rl = AdvancedRootLocusPlotter(num, den) roots = rl.compute_root_locus() # 计算关键性能指标 damping_ratios = [] natural_freqs = [] for root_set in roots: for root in root_set: if np.imag(root) != 0: # 只考虑复数极点 wn = np.abs(root) zeta = -np.real(root) / wn damping_ratios.append(zeta) natural_freqs.append(wn) # 可视化性能参数 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6)) # 阻尼比变化 ax1.plot(rl.k_values[:len(damping_ratios)], damping_ratios, 'g-') ax1.set_xscale('log') ax1.set_xlabel('K值') ax1.set_ylabel('阻尼比ζ') ax1.grid(True) # 自然频率变化 ax2.plot(rl.k_values[:len(natural_freqs)], natural_freqs, 'r-') ax2.set_xscale('log') ax2.set_xlabel('K值') ax2.set_ylabel('自然频率ω_n') ax2.grid(True) plt.tight_layout() plt.show() # 输出设计建议 optimal_idx = np.argmin(np.abs(np.array(damping_ratios) - 0.707)) # 最佳阻尼比 print(f"推荐K值范围: {rl.k_values[optimal_idx]:.2f} 附近") print(f"对应阻尼比: {damping_ratios[optimal_idx]:.3f}") print(f"自然频率: {natural_freqs[optimal_idx]:.2f} rad/s") # 分析示例系统 system_performance_analysis(num=[1, 1], den=[1, 5, 6, 0])

这个工具能帮助我们：

观察不同K值下的阻尼比变化
确定最佳响应区域（通常ζ≈0.707）
避免进入不稳定区域（极点进入右半平面）

6. 工程扩展：处理复杂系统与性能优化

实际工程系统往往更复杂。我们最后实现一个能处理以下情况的增强版本：

非最小相位系统（右半平面零点）
多参数调节（同时变化多个极点/零点位置）
灵敏度分析（参数微小变化的影响）

class EngineeringRootLocusPlotter(AdvancedRootLocusPlotter): def __init__(self, num, den, variable_param=None): """ variable_param: 可变参数配置，如{'position': 'zero', 'index': 0, 'range': [1, 10]} """ super().__init__(num, den) self.variable_param = variable_param def compute_variable_locus(self): """计算参数变化时的根轨迹簇""" if not self.variable_param: return self.compute_root_locus() param_values = np.linspace(*self.variable_param['range'], 10) all_roots = [] for val in param_values: # 动态修改系统参数 modified_num, modified_den = self._modify_system(val) # 计算当前参数下的根轨迹 temp_roots = [] for k in self.k_values: closed_loop_den = modified_den + k * np.pad(modified_num, (len(modified_den)-len(modified_num), 0), 'constant') root = np.roots(closed_loop_den) temp_roots.append(root) all_roots.append(np.array(temp_roots)) return np.array(all_roots) def _modify_system(self, param_value): """根据参数值修改系统零极点""" num = self.num.copy() den = self.den.copy() if self.variable_param['position'] == 'zero': roots = np.roots(num) roots[self.variable_param['index']] = -param_value num = np.poly(roots) elif self.variable_param['position'] == 'pole': roots = np.roots(den) roots[self.variable_param['index']] = -param_value den = np.poly(roots) return num, den # 使用示例：观察零点位置变化的影响 eng_rl = EngineeringRootLocusPlotter( num=[1, 2], den=[1, 5, 6, 0], variable_param={'position': 'zero', 'index': 0, 'range': [1, 5]} ) root_clusters = eng_rl.compute_variable_locus() # 可视化参数变化影响 fig, ax = plt.subplots(figsize=(10, 8)) colors = plt.cm.viridis(np.linspace(0, 1, len(root_clusters))) for i, roots in enumerate(root_clusters): for j in range(roots.shape[1]): ax.plot(roots[:, j].real, roots[:, j].imag, color=colors[i], alpha=0.7) ax.set_xlabel('Real Axis') ax.set_ylabel('Imaginary Axis') ax.grid(True) plt.show()

这种分析特别适用于：