在对动态系统进行评估时会经常使用阶跃响应,这是因为阶跃输入对系统来说是最严峻的工作状态,如果在阶跃函数作用下的动态性能满足要求,那么系统在其它形式的函数作用下,其动态性能也能够得到保证。
线性时不变 (LTI-linear time invariant) 系统可以等效地描述为传递函数、状态空间模型,或使用 ODE 积分器进行数值求解,之前已经对这三类方法分别进行了介绍。今天针对不同系统的阶跃响应来回顾这三类方法在Python中的实现。
一阶系统
假设一个一阶线性微分方程:
其中,tp、Kp为常数;u为输入;y为输出。
这里可以使用传递函数、 状态空间模型和半显式微分方程的方法表示此微分方程。
传递函数
状态空间模型
微分方程
这里假设Kp=5,tp=3。
Python代码:
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from scipy.integrate import odeint
Kp = 5.0
tp = 3.0
# 传递函数
num = [Kp]
den = [tp,1]
sys1 = signal.TransferFunction(num,den)
t1,y1 = signal.step(sys1)
# 状态空间
A = -1.0/tp
B = Kp/tp
C = 1.0
D = 0.0
sys2 = signal.StateSpace(A,B,C,D)
t2,y2 = signal.step(sys2)
# 微分方程
def model3(y,t):
u = 1
return (-y + Kp * u)/tp
t3 = 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()
二阶系统
对于典型的二阶系统,响应取决于它是过阻尼、临界阻尼还是欠阻尼二阶系统。
微分方程
其中,Kp为增益,ζ为阻尼系数,τ为二阶系统时间常数,θp为死区时间。
传递函数
状态空间
这里令Kp=3 , τ=1 ,ζ=0.25,θ=0
Python代码:
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from scipy.integrate import odeint
Kp = 3.0 # gain
tau = 1.0 # time constant
zeta = 0.25 # damping factor
theta = 0.0 # no time delay
du = 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.0
sys2 = 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 * du
plt.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 ——
联系客服