用PythonSymPy玩转自动控制原理从微分方程到传递函数的代码实践刚接触自动控制原理时你是否曾被那些复杂的微分方程和拉氏变换搞得头晕眼花传统教材往往聚焦于理论推导却忽略了如何将这些数学工具真正用起来。今天我将带你用Python的SymPy库把课本上的抽象公式转化为可执行的代码让数学建模变得触手可及。1. 环境准备与SymPy基础工欲善其事必先利其器。我们需要一个能处理符号计算的工具——SymPy。这个纯Python库无需额外依赖安装简单pip install sympySymPy的核心功能包括符号变量定义像数学一样声明未知数方程求解代数方程、微分方程都能解拉氏变换内置完整的积分变换工具链表达式简化自动化简复杂数学表达式先来感受下基本操作from sympy import symbols, Eq, Function, dsolve t symbols(t) # 定义时间变量 y Function(y)(t) # 定义未知函数y(t) dy y.diff(t) # y对t的一阶导数 equation Eq(dy 3*y, 2) # 创建微分方程dy/dt 3y 2 solution dsolve(equation) # 求解微分方程 print(solution) # 输出解析解这段代码求解了一个简单的一阶线性微分方程SymPy会返回包含积分常数的通解。相比手工计算我们节省了推导时间还能确保结果绝对准确。2. 微分方程建模实战控制系统的核心是微分方程模型。让我们通过一个经典案例——弹簧-质量-阻尼系统看看如何用代码实现全过程。2.1 系统建模考虑一个质量为m的物体连接弹簧(k)和阻尼器(c)受力F(t)作用from sympy import symbols, Function, Eq, dsolve t, m, c, k symbols(t m c k) # 定义符号变量 x Function(x)(t) # 位移函数 F Function(F)(t) # 外力函数 # 牛顿第二定律建立微分方程 eq Eq(m*x.diff(t,2) c*x.diff(t) k*x, F) print(系统微分方程:, eq)运行后会输出标准的二阶微分方程形式m⋅d²x/dt² c⋅dx/dt k⋅x F(t)2.2 方程求解假设参数m1, c2, k1外力F(t)sin(t)初始条件x(0)0, x(0)1from sympy import sin # 代入具体参数 eq_sub eq.subs({m:1, c:2, k:1, F:sin(t)}) # 求解带初始条件的微分方程 solution dsolve(eq_sub, ics{x.subs(t,0):0, x.diff(t).subs(t,0):1}) print(解析解:, solution.simplify())SymPy会返回包含三角函数和指数项的完整解析解。如果想看数值结果可以进一步用lambdify将符号解转为可计算函数from sympy import lambdify import numpy as np import matplotlib.pyplot as plt x_func lambdify(t, solution.rhs, numpy) t_vals np.linspace(0, 10, 100) x_vals x_func(t_vals) plt.plot(t_vals, x_vals) plt.xlabel(Time (s)) plt.ylabel(Displacement (m)) plt.title(System Response) plt.grid(True) plt.show()3. 拉氏变换与传递函数拉氏变换将时域微分方程转换为频域代数方程极大简化了系统分析。SymPy的laplace_transform模块让这些变换变得轻而易举。3.1 基本变换操作计算常见函数的拉氏变换from sympy import exp, Heaviside, laplace_transform s symbols(s) a, t symbols(a t, positiveTrue) # 单位阶跃函数 print(laplace_transform(Heaviside(t), t, s)[0]) # 指数函数 print(laplace_transform(exp(-a*t), t, s)[0]) # 正弦函数 from sympy import sin print(laplace_transform(sin(a*t), t, s)[0])3.2 微分方程的拉氏变换解法重新考虑前面的弹簧-质量-阻尼系统这次用拉氏变换法求解from sympy import laplace_transform, inverse_laplace_transform # 对微分方程两边进行拉氏变换 L_x symbols(L_x) transformed_eq eq_sub.replace( lambda expr: expr.is_Derivative, lambda expr: laplace_transform(expr, t, s)[0] ).replace(x, L_x).replace(laplace_transform(sin(t), t, s)[0], 1/(s**21)) # 解代数方程求L_x L_x_sol solve(transformed_eq, L_x)[0] # 拉氏反变换得到时域解 x_sol inverse_laplace_transform(L_x_sol, s, t) print(拉氏变换法得到的解:, x_sol.simplify())3.3 传递函数推导传递函数是零初始条件下输出与输入的拉氏变换比from sympy import solve # 假设初始条件为零 transformed_eq_zero_ics transformed_eq.subs({x.subs(t,0):0, x.diff(t).subs(t,0):0}) G solve(transformed_eq_zero_ics, L_x)[0] / (1/(s**21)) # 输出/输入 print(传递函数:, G.simplify())这会输出标准的二阶系统传递函数形式1/(s**2 2*s 1)4. 系统分析与可视化有了传递函数我们可以进一步分析系统特性。4.1 零极点分析from sympy import roots, plot_poles # 提取分母多项式得到特征方程 denominator s**2 2*s 1 # 计算极点 poles roots(denominator, s) print(系统极点:, poles) # 绘制零极点图 from sympy.plotting import plot_pole_zero G 1/denominator plot_pole_zero(G, s, showTrue)4.2 频率响应分析虽然SymPy不直接支持Bode图但我们可以结合scipy和matplotlib实现from scipy import signal import matplotlib.pyplot as plt # 传递函数系数 (s² 2s 1) num [1] den [1, 2, 1] # 创建系统 sys signal.TransferFunction(num, den) # 绘制Bode图 w, mag, phase signal.bode(sys) plt.figure() plt.semilogx(w, mag) # Bode magnitude plot plt.title(Bode Plot - Magnitude) plt.figure() plt.semilogx(w, phase) # Bode phase plot plt.title(Bode Plot - Phase) plt.show()4.3 时域响应模拟t, y signal.step(sys) plt.plot(t, y) plt.title(Step Response) plt.xlabel(Time [s]) plt.ylabel(Amplitude) plt.grid(True) plt.show()5. 进阶应用状态空间分析对于更复杂的系统状态空间表示更为强大。SymPy也能处理from sympy import Matrix # 状态空间矩阵 A Matrix([[0, 1], [-1, -2]]) B Matrix([[0], [1]]) C Matrix([[1, 0]]) D Matrix([0]) # 计算状态转移矩阵 from sympy import eye s symbols(s) Phi (s*eye(2) - A).inv() print(状态转移矩阵:, Phi) # 传递函数矩阵 G C * Phi * B D print(传递函数:, G[0].simplify())6. 实际工程问题解决最后我们看一个实际的电机控制系统例子。直流电机的微分方程为# 电机参数 J, b, K, R, L symbols(J b K R L) # 状态变量电流i和角速度ω i Function(i)(t) ω Function(ω)(t) V Function(V)(t) # 电学方程 eq1 Eq(L*i.diff(t) R*i K*ω, V) # 力学方程 eq2 Eq(J*ω.diff(t) b*ω, K*i) # 求解传递函数ω/V from sympy import solve sols solve((eq1, eq2), (i.diff(t), ω.diff(t)))通过类似前面的方法可以推导出电机转速与输入电压之间的传递函数进而分析系统性能。