【原】基于Python的动态建模——系统的阶跃响应

在对动态系统进行评估时会经常使用阶跃响应，这是因为阶跃输入对系统来说是最严峻的工作状态，如果在阶跃函数作用下的动态性能满足要求，那么系统在其它形式的函数作用下，其动态性能也能够得到保证。

线性时不变 (LTI-linear time invariant) 系统可以等效地描述为传递函数、状态空间模型，或使用 ODE 积分器进行数值求解，之前已经对这三类方法分别进行了介绍。今天针对不同系统的阶跃响应来回顾这三类方法在Python中的实现。

01

一阶系统

假设一个一阶线性微分方程：

其中，tp、Kp为常数；u为输入；y为输出。

这里可以使用传递函数、 状态空间模型和半显式微分方程的方法表示此微分方程。

• 传递函数

• 状态空间模型
  

• 微分方程
  

这里假设Kp=5，tp=3。

Python代码：

import numpy as npfrom scipy import signalimport matplotlib.pyplot as pltfrom scipy.integrate import odeint
Kp = 5.0tp = 3.0
# 传递函数num = [Kp]den = [tp,1]sys1 = signal.TransferFunction(num,den)t1,y1 = signal.step(sys1)
# 状态空间A = -1.0/tpB = Kp/tpC = 1.0D = 0.0sys2 = signal.StateSpace(A,B,C,D)t2,y2 = signal.step(sys2)
# 微分方程def model3(y,t):    u = 1    return (-y + Kp * u)/tpt3 = np.linspace(0,14,100)y3 = odeint(model3,0,t3)
plt.figure(1)plt.plot(t1,y1,'r--',linewidth=3,label='Transfer Function')plt.plot(t2,y2,'g:',linewidth=2,label='State Space model')plt.plot(t3,y3,'b-',linewidth=1,label='ODE Integrator')plt.xlabel('Time')plt.ylabel('Response (y)')plt.legend(loc='best')plt.show()

02

二阶系统

对于典型的二阶系统，响应取决于它是过阻尼、临界阻尼还是欠阻尼二阶系统。

• 微分方程
  

其中，Kp为增益，ζ为阻尼系数，τ为二阶系统时间常数，θp为死区时间。

• 传递函数

• 状态空间
  

这里令Kp=3 ， τ=1 ，ζ=0.25，θ=0

Python代码：

import numpy as npfrom scipy import signalimport matplotlib.pyplot as pltfrom scipy.integrate import odeint
Kp = 3.0    # gaintau = 1.0   # time constantzeta = 0.25 # damping factortheta = 0.0 # no time delaydu = 1.0    # change in u
# 传递函数num = [Kp]den = [tau**2,2*zeta*tau,1]sys1 = signal.TransferFunction(num,den)t1,y1 = signal.step(sys1)
# 状态空间A = [[0.0,1.0],[-1.0/tau**2,-2.0*zeta/tau]]B = [[0.0],[Kp/tau**2]]C = [1.0,0.0]D = 0.0sys2 = signal.StateSpace(A,B,C,D)t2,y2 = signal.step(sys2)
# 微分方程def model3(x,t):    y = x[0]    dydt = x[1]    dy2dt2 = (-2.0*zeta*tau*dydt - y + Kp*du)/tau**2    return [dydt,dy2dt2]t3 = np.linspace(0,25,100)x3 = odeint(model3,[0,0],t3)y3 = x3[:,0]
plt.figure(1)plt.plot(t1,y1*du,'r--',linewidth=3,label='Transfer Fcn')plt.plot(t2,y2*du,'g:',linewidth=2,label='State Space')plt.plot(t3,y3,'b-',linewidth=1,label='ODE Integrator')y_ss = Kp * duplt.plot([0,max(t1)],[y_ss,y_ss],'k:')plt.xlim([0,max(t1)])plt.xlabel('Time')plt.ylabel('Response (y)')plt.legend(loc='best')plt.show()

—— end ——