基于Python的动态建模——传递函数的实现与应用

传递函数是动态系统的输入到输出的表示。频域（相对于时域）的一个好处是将微分方程转换为代数方程。这些代数方程可以重新排列并转换回时域以获得解决方案，或者进一步与其他传递函数相结合以创建更复杂的系统。创建传递函数的第一步是用拉普拉斯变换来转换微分方程。

传递函数G(s)将输入U(s)与输出Y(s)联系起来：

01

使用Python进行拉普拉斯变换

要得到传递函数，首先要对系统的输出与输入量进行拉普拉斯变换。

这里要使用到Sympy，Sympy是一个支持符号计算的Python库，符号计算是指处理数学对象的计算，与数值计算相比，在符号计算中，数学对象是精确表示的，而不是近似的，未计算的数学表达式会以符号形式保留。

典型的符号计算有：

• 表达式化简
• 表达式求值
• 表达式的变形：展开、积、幂、部分分式表示、将三角函数转换为指数函数等
• 一元或多元微分
• 带条件的化简
• 部分或完整的因式分解
• 求解线性或非线性方程
• 求解微分方程或差分方程
• 求极限
• 求函数的定积分、不定积分
• 泰勒展开、洛朗展开等
• 无穷级数展开
• 级数求和
• 矩阵运算

Sympy同样支持拉普拉斯变换(laplace transform)、拉普拉斯逆变换(inverse laplace transform)、部分分数展开( apart)、多项式展开(expand)和多项式根(roots)等系统建模分析中常用的运算。

下面以常见的拉普拉斯变换为例进行介绍：

import sympy as symfrom sympy.abc import s,t,x,y,z,afrom sympy.integrals import laplace_transformfrom sympy.integrals import inverse_laplace_transform

• 单位阶跃

U = laplace_transform(1,t,s)U[0]

• 

• sin(at)
  

02

计算案例

这里假设一包含阶跃和斜坡函数的输入，其图像为：

系统传递函数是一个具有两个极点（分母根）和一个零点（分子根）的稳定系统：

• 基于Sympy的符号解

首先将阶跃和斜坡信号分解为三个单独的函数。分别计算系统对这三个输入的响应，然后对信号求和。

这里利用到了延迟定理（t域平移定理）：

求解代码：

import sympy as symfrom sympy.abc import s,t,x,y,zimport numpy as npfrom sympy.integrals import inverse_laplace_transformimport matplotlib.pyplot as plt
# 定义输入# 阶跃U1 = 2/s*sym.exp(-s)# 3s处的斜坡U2 = -1/s**2*sym.exp(-3*s)# 5s处斜坡U3 = 1/s**2*sym.exp(-5*s)
# 定义传递函数G = 5*(s+1)/(s+3)**2
# 计算输出Y1 = G * U1Y2 = G * U2Y3 = G * U3
# 拉普拉斯逆变换u1 = inverse_laplace_transform(U1,s,t)u2 = inverse_laplace_transform(U2,s,t)u3 = inverse_laplace_transform(U3,s,t)y1 = inverse_laplace_transform(Y1,s,t)y2 = inverse_laplace_transform(Y2,s,t)y3 = inverse_laplace_transform(Y3,s,t)print('y1')print(y1)
# 创建绘图数据tm = np.linspace(0,8,100)us = np.zeros(len(tm))ys = np.zeros(len(tm))
# 数值化u和yfor u in [u1,u2,u3]:    for i in range(len(tm)):        us[i] += u.subs(t,tm[i])for y in [y1,y2,y3]:    for i in range(len(tm)):        ys[i] += y.subs(t,tm[i])
# 结果图plt.figure()plt.plot(tm,us,label='u(t)')plt.plot(tm,ys,'--r',label='y(t)')plt.legend()plt.xlabel('Time')plt.show()

• 基于GEKKO的数值解法
  

符号解法的一个替代方法是在时域中用数值法计算响应。传递函数必须首先被转化为微分方程。

from gekko import GEKKOimport numpy as npimport matplotlib.pyplot as plt
m = GEKKO()
# 创建数据点nt = 81m.time = np.linspace(0,8,nt)
# 定义输入ut = np.zeros(nt)# 阶跃ut[11:31] = 2.0# 斜坡for i in range(31,51):    ut[i] = ut[i-1] - 0.1
# 定义模型u = m.Param(value=ut)ud = m.Var()y = m.Var()dydt = m.Var()m.Equation(ud==u)m.Equation(dydt==y.dt())m.Equation(dydt.dt() + 6*y.dt() + 9*y==5*ud.dt()+5*u)
# 计算设置m.options.IMODE=7m.options.NODES=4m.solve(disp=False)
#结果图plt.figure()plt.plot(m.time,u.value,label='u(t)')plt.plot(m.time,y.value,'--g',label='y(t)')plt.legend()plt.xlabel('Time')plt.show()

符号（解析）解的优点是高度准确，不依赖数值方法来近似解。此外，解决方案是一个紧凑的形式，可以用于进一步的分析。符号解仅限于输入函数和系统传递函数可以用拉普拉斯变换的情况。拉普拉斯变换的符号解法也无法用于非线性或复杂的系统

数值解可以处理成千上万或数百万的非线性关系的方程。数值解法的缺点是它是对真实解法的一种近似，可能有不准确的地方。另一个缺点是，求解器可能会无法收敛。