Songyuli's blog
0%

被控对象的行为

Posted on Edited on In Disqus:

被控对象的行为

准备工作

为了更加方便地绘制相关图像,在这里要进行相应的准备工作

1
2
3
4
import matplotlib.pyplot as plt
import numpy as np
import control as ct 
from control.matlab import *

绘图时用以确定线条类型

1
2
3
4
5
6
def line_generator():
    linestyle=["-","--","-.",":"]
    lineid=0
    while True:
        yield linestyle[lineid] #* 这个yield用法很有趣,可以学习一下
        lineid=(lineid+1)%len(linestyle)

用以完善绘图

1
2
3
4
5
6
def plot_set(fig_ax,*args):
    fig_ax.set_xlabel(args[0])
    fig_ax.set_ylabel(args[1])
    fig_ax.grid(ls=':')
    if len(args)==3:
        fig_ax.legend(loc=args[2])

用以完善伯德图函数

1
2
3
4
5
6
7
8
9
10
def bodeplot_set(fig_ax,*args):
    fig_ax[0].grid(which='both',ls=':')
    fig_ax[0].set_ylabel('Gain [dB]')
    fig_ax[1].grid(which='both',ls=':')
    fig_ax[1].set_xlabel('$\omega$ [rad/s]')
    fig_ax[1].set_ylabel('Phase [deg]')
    if len(args)>0:
        fig_ax[1].legend(loc=args[0])
    if len(args)>1:
        fig_ax[0].legend(loc=args[1])

时域响应

在这一部分,我们选取常见的阶跃函数step对于一阶和二阶的滞后系统进行描述

一阶滞后系统

之前用来描述有阻尼手推车系统的传递函数的传递函数为
$$\mathcal{P}(s)=\frac{\frac{1}{\mu}}{1+\frac{M}{\mu}s}$$
可以将$\frac{1}{\mu}$简化为增益$K$,$\frac{M}{\mu}$简化为时间常数$T$,
这一类系统的阶跃响应可以被表示为

1
2
3
4
5
6
7
ts,k=0.5,1
p=tf([0,k],[ts,1])
y,t=step(p,np.linspace(0,5,100))

fig,ax=plt.subplots()
ax.plot(t,y)
plot_set(ax,'t','y')

​
png
​

其中,我们需要注意到,$T$代表系统达到稳定值的63.2%所需要的时间.
下面展示了在不同的$T$下系统的变化

1
2
3
4
5
6
7
8
9
10
ls=line_generator()
fig,ax=plt.subplots()

k=1
tl=[0.1,0.5,1]
for tt in tl:
    p=tf([0,k],[tt,1])
    y,t=step(p,np.linspace(0,5,100))
    ax.plot(t,y,next(ls),label="T="+str(tt))
plot_set(ax,"t","y","best")

​
png
​

可以明显看出$T$越大,达到平衡就越慢,$K$为增益,在这里不再加以演示了

二阶滞后系统

如果考量一个阻尼弹簧的模型,我们可以得到其传递模型为
$$\mathcal{P}(s)=\frac{\omega_n^2}{s^2+2\zeta\omega_ns+\omega_n^2}$$
式中$\zeta$被称为阻尼系数,$\omega_n$被称为无阻尼自然频率.
这样的系统被称为二阶滞后系统

实际上就是阻尼振动(LCR振荡的模型)

1
2
3
4
5
6
7
zeta,omega0=0.4,5
p=tf([0,omega0**2],[1,2*zeta*omega0,omega0**2])
y,t=step(p,np.linspace(0,5,100))

fig,ax=plt.subplots()
ax.plot(t,y)
plot_set(ax,"t","y")

​
png
​

可以注意到响应图象呈现先增大后减小(中间存在过冲)过冲在二阶滞后系统时有发生.
下面展现出过冲在不同的阻尼系数下的表现

1
2
3
4
5
6
7
8
9
10
ls=line_generator()
fig,ax=plt.subplots()

omega=5
beta=[0.3,0.5,1]
for tt in beta:
    p=tf([0,omega**2],[1,2*tt*omega,omega**2])
    y,t=step(p,np.linspace(0,5,100))
    ax.plot(t,y,next(ls),label="$\zeta$="+str(tt))
plot_set(ax,"t","y","best")

​
png
​

可以注意到随着$\zeta$增加,过冲逐渐减小,乃至消失,而不同的$\omega_n$的影响可以如下图所示

1
2
3
4
5
6
7
8
9
10
ls=line_generator()
fig,ax=plt.subplots()

omegas=[1,5,10]
tt=0.3
for omega in omegas:
    p=tf([0,omega**2],[1,2*tt*omega,omega**2])
    y,t=step(p,np.linspace(0,5,100))
    ax.plot(t,y,next(ls),label="$\omega_n$="+str(omega))
plot_set(ax,"t","y","best")

​
png
​

可见$\omega_n$不会改变响应的大小,只会改变响应速度

下面展现了一个高阶响应系统的阶跃响应

1
2
3
4
5
6
p=tf([0,1],[1,2,3,1])
y,t=step(p,np.linspace(0,5,100))

fig,ax=plt.subplots()
ax.plot(t,y)
plot_set(ax,"t","y")

​
png
​

状态空间模型的时域响应

在状态空间模型中,我们还会考虑不同的初始值会给出什么样不同的结果

状态空间模型(单输入,单输出)定义
$$\mathbf{\dot{x}}(t)=\mathrm{A}\mathbf{x}(t)+\mathbf{B}u(t)$$
$$y(t)=\mathrm{C}\mathbf{x}(t)+Du(t)$$
我们考虑一个系统,其中$A=\left[\begin{matrix}0&1\-4&-5\end{matrix}\right],B=\left[\begin{matrix}0\1\end{matrix}\right],C=\left[\begin{matrix}1&0\0&1\end{matrix}\right],D=\left[\begin{matrix}0\0\end{matrix}\right]$,为了便于观察中间状态的行为,将$C$设定为单位阵
而在Python中,可以利用initial函数获得系统对于初值的响应

1
2
3
4
5
6
7
8
9
10
11
12
A=[[0,1],[-4,-0.5]]
B=[[0],[1]]
C=[[1,0],[0,1]]
D=[[0],[0]]
P=ss(A,B,C,D)
x0=[-0.3,0.4]
x,t=initial(P,np.linspace(0,30,1000),x0)

fig,ax=plt.subplots()
ax.plot(t,x[:,0],label='$x_1$')
ax.plot(t,x[:,1],'-.',label='$x_2$')
plot_set(ax,"t","x",'best')

​
png
​

可以注意到,两个状态(分别对应于阻尼运动的位置和速度)最终都会趋向0,这是阻尼运动的特征
对于一个微分方程$\dot{x}(t)=ax(t)$,其具有解$x(t)=e^{at}x(0)$,而将$a$替换为矩阵$\mathrm{A}$可以得到状态方程的解为
$$\mathbf{x}(t)=e^{\mathrm{A}t}\mathbf{x}(0)$$
如果可以计算出这个指数,那么就可以得到状态方程的行为特征

1
2
3
4
import scipy.linalg
mA=np.array(A)
t=10
scipy.linalg.expm(mA*t)
array([[ 0.05345953,  0.03466486],
       [-0.13865943,  0.0361271 ]])

接下来我们考虑具有输入的情况,利用下面的求解公式
$$\mathrm{x}(t)=e^{\mathrm{A}t}\mathrm{x}(0)+\int_0^te^{\mathrm{A}(t-t_s)}\mathrm{B}u(t_s)\mathrm{d}t_s$$
其中右边第一项为零输入响应,第二项为零初值响应
接下来来看一下之前函数$x=\left[\begin{matrix}0\0\end{matrix}\right]$时的零初值响应,输入为阶跃输入

1
2
3
4
5
6
7
Td=np.linspace(0,20,1000)
x,t=step(P,Td)

fig,ax=plt.subplots()
ax.plot(t,x[:,0],label="$x_1$")
ax.plot(t,x[:,1],ls="-.",label="$x_2$")
plot_set(ax,"t","x",'best')

​
png
​

我们可以进一步将之前的零输入响应添加进去,只要将两者相加即可,或者使用lsim

1
2
3
4
5
6
7
8
9
Td=np.linspace(0,20,1000)
x1,t=step(P,Td)
x2,t=initial(P,Td,[0.3,-0.4])
x=x1[...,0]+x2

fig,ax=plt.subplots()
ax.plot(t,x[:,0],label="$x_1$")
ax.plot(t,x[:,1],ls="-.",label="$x_2$")
plot_set(ax,"t","x",'best')

​
png
​

下面我们可以稍稍修改一下代码,实现$u(t)=3\sin 5t,x(0)=[0.5,1]^T$时系统响应

对于一般的响应,建议还是使用lsim实现

1
2
3
4
5
6
7
8
9
Td=np.linspace(0,20,1000)
x1=3*np.sin(5*Td)
x2,t=initial(P,Td,[0.5,1])
x,_,_=lsim(P,x1,Td,[0.5,1])

fig,ax=plt.subplots()
ax.plot(t,x[:,0],label="$x_1$")
ax.plot(t,x[:,1],ls="-.",label="$x_2$")
plot_set(ax,"t","x",'best')

​
png
​

稳定性

在研究系统工作时,我们发现一部分值可能会导致系统发散,即系统不稳定,我们需要研究系统的稳定性取决于哪些因素

输入输出稳定性

如果输入信号有界,输出也有界,这种情况被称作输入输出稳定(BIBO稳定)

判断不稳定的方法
观察传递函数$\mathcal{P}(s)$的极点(其中代数式达到无穷的情况)
在python中可以使用pole函数寻找到传递函数的极点

1
2
P1=tf([0,1],[1,1])
print("P1:",pole(P1))
P1: [-1.+0.j]

接下来看看二阶滞后函数中$K=1,\omega_n=1,\zeta=0.025$时的传递函数的极点

1
2
P2=tf([0,1],[1,0.05,1])
print("P2:",pole(P2))
P2: [-0.025+0.99968745j -0.025-0.99968745j]

系统处于输入输出稳定条件为:传递函数的所有极点实部均为负数
实部为负值的极点为稳定极点,否则为不稳定极点

渐进稳定性

前面介绍的是相对于传递函数而言的,现在讨论状态空间模型,可以通过观察系统的A矩阵的特征值实现对稳定性的观测,则有:
系统稳定的充分必要条件:矩阵A所有特征值的实部均为负数

但是同时需要注意到的是这里的稳定性具有渐进的特性即
$$\lim_{t\to\infin}x(t)=0$$
如果具有渐进稳定性,那么对于有界输入,存在有界输出(逆命题不成立)

使用np.linalg.eigvals可以求得矩阵的特征值

1
2
A=np.array([[0,1],[-4,-5]])
np.linalg.eigvals(A)
array([-1., -4.])

而通过绘制相平面图(即状态$x$的演化轨迹)如下代码所示

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def make_phase(matA,w=1.5):
    Y,X=np.mgrid[-w:w:100j,-w:w:100j]
    s,v=np.linalg.eig(matA)
    print(s)

    U=matA[0,0]*X+matA[0,1]*Y
    V=matA[1,0]*X+matA[1,1]*Y

    t=np.arange(-k,k,0.01)
    fig,ax=plt.subplots()
    if s.imag[0]==0 and s.imag[1]==0:
        ax.plot(t,(v[1,0]/v[0,0])*t,ls="-")
        ax.plot(t,(v[1,1]/v[0,1])*t,ls="-")
    ax.streamplot(X,Y,U,V,density=0.7,color='k')
    plot_set(ax,'$x_1$','$x_2')

先来刚刚那个矩阵

1
make_phase(A,w=4)
[-1. -4.]

png

可以注意到,所有的箭头都指向不动线的,下面展示如果特征值大于零会如何

1
make_phase(np.array([[0,1],[1,4]]),w=4)
[-0.23606798  4.23606798]

png

还有没有实数特征值的

1
make_phase(np.array([[np.cos(0.5),np.sin(0.5)],[-np.sin(0.5),np.cos(0.5)]]),w=4)
[0.87758256+0.47942554j 0.87758256-0.47942554j]

png

这些极点和系统行为密切相关,若极点负值越大,那么响应越迅速

当极点的虚部不为零时,会出现振荡,虚部越大振荡越快

观察实部,可以得到振荡振幅;观察虚部,可以得到振荡周期

频域响应

当输入为冲激输入(函数形式为$\delta$函数)时,系统的响应就是传递函数,因此只需要对施加冲激输入后的系统响应进行拉普拉斯变换就可以得到传递函数

但是在现实状况下,实现冲激输入非常困难,因此我们需要对信号进行傅里叶变换,研究输入输出的频域响应得到结果

我们可以观察二阶滞后系统在输入正弦波时的响应

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
fig,ax=plt.subplots(2,2)

zeta=0.7
omega_n=5
P=tf([0,omega_n**2],[1,2*zeta*omega_n,omega_n**2])

freq=[2,5,10,20]
Td=np.linspace(0,5,1000)
for i in range(2):
    for j in range(2):
        u=np.sin(freq[2*i+j]*Td)
        y,t,x0=lsim(P,u,Td,0)
        ax[i,j].plot(t,u,ls='--',label='u')
        ax[i,j].plot(t,y,label='y')
        plot_set(ax[i,j],'t','u,y')
ax[0,0].legend()
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\timeresp.py:935: UserWarning: return_x specified for a transfer function system. Internal conversion to state space used; results may meaningless.
  warnings.warn(
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\timeresp.py:935: UserWarning: return_x specified for a transfer function system. Internal conversion to state space used; results may meaningless.
  warnings.warn(
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\timeresp.py:935: UserWarning: return_x specified for a transfer function system. Internal conversion to state space used; results may meaningless.
  warnings.warn(
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\timeresp.py:935: UserWarning: return_x specified for a transfer function system. Internal conversion to state space used; results may meaningless.
  warnings.warn(





<matplotlib.legend.Legend at 0x1aa4e6be160>

​
png
​

可以注意到随着频率逐渐增大,振幅逐渐减小,同时相位的滞后就会比较大

对于各个频率,振幅比以$20\log_10\frac{B(\omega)}{A}$(分贝,dB)表示,绘制出图象为幅频图
将相位(deg)绘出的图形称为相频图,两者合称伯德图

可以使用bode获得伯德图(记得转换)

一阶滞后系统的频域响应

下面我们来研究一阶滞后系统的伯德图,我们分别取时间常数$T=1,0.5,0.1$

1
2
3
4
5
6
7
8
9
10
11
K=1
T=[1,0.5,0.1]
Ls=line_generator()
fig,ax=plt.subplots(2,1)
for i in range(len(T)):
    P=tf([0,K],[T[i],1])
    gain,phase,w=bode(P,logspace(-2,2),Plot=False)
    pltargs={'ls':next(Ls),'label':'T='+str(T[i])}
    ax[0].semilogx(w,20*np.log10(gain),**pltargs)
    ax[1].semilogx(w,phase*180/np.pi,**pltargs)
bodeplot_set(ax,3,3)
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\freqplot.py:187: FutureWarning: 'Plot' keyword is deprecated in bode_plot; use 'plot'
  warnings.warn("'Plot' keyword is deprecated in bode_plot; use 'plot'",
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\freqplot.py:187: FutureWarning: 'Plot' keyword is deprecated in bode_plot; use 'plot'
  warnings.warn("'Plot' keyword is deprecated in bode_plot; use 'plot'",
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\freqplot.py:187: FutureWarning: 'Plot' keyword is deprecated in bode_plot; use 'plot'
  warnings.warn("'Plot' keyword is deprecated in bode_plot; use 'plot'",

png

可以注意到,随着频率的提高,输出信号会发生衰减和相移,时间常数$T$越大,衰减和相移越为明显

这很好理解,时间常数越大,代表响应越慢,自然衰减和相移就大

一般可以认为,在$\frac{1}{T}$的范围之内,输入信号振幅与输出信号相似(此时相移45°)

二阶滞后系统的频域响应

在这一部分,我们来研究一下二阶滞后系统的频域响应,取阻尼系数$\zeta=1,0.7,0.4$

1
2
3
4
5
6
7
8
9
10
11
zeta=[1,0.7,0.4]
omega=1
Ls=line_generator()
fig,ax=plt.subplots(2,1)
for i in range(len(T)):
    P=tf([0,omega**2],[1,2*zeta[i]*omega,omega**2])
    gain,phase,w=bode(P,logspace(-2,2),Plot=False)
    pltargs={'ls':next(Ls),'label':'$\zeta$='+str(T[i])}
    ax[0].semilogx(w,20*np.log10(gain),**pltargs)
    ax[1].semilogx(w,phase*180/np.pi,**pltargs)
bodeplot_set(ax,3,3)
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\freqplot.py:187: FutureWarning: 'Plot' keyword is deprecated in bode_plot; use 'plot'
  warnings.warn("'Plot' keyword is deprecated in bode_plot; use 'plot'",
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\freqplot.py:187: FutureWarning: 'Plot' keyword is deprecated in bode_plot; use 'plot'
  warnings.warn("'Plot' keyword is deprecated in bode_plot; use 'plot'",
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\freqplot.py:187: FutureWarning: 'Plot' keyword is deprecated in bode_plot; use 'plot'
  warnings.warn("'Plot' keyword is deprecated in bode_plot; use 'plot'",

png

通过观察伯德图,我们可以得到在高频段振幅会逐渐减小,虽然在较高频下振幅的衰减与$\zeta$无关,但是如果$\zeta$较低,会出现过冲(类似于共振)
相位在二姐滞后系统最终会移相180°

我们再来看看改变频率$\omega_n$的结果

1
2
3
4
5
6
7
8
9
10
11
zeta=0.5
omega=[1,5,10]
Ls=line_generator()
fig,ax=plt.subplots(2,1)
for i in range(len(T)):
    P=tf([0,omega[i]**2],[1,2*zeta*omega[i],omega[i]**2])
    gain,phase,w=bode(P,logspace(-2,2),Plot=False)
    pltargs={'ls':next(Ls),'label':'$\omega$='+str(T[i])}
    ax[0].semilogx(w,20*np.log10(gain),**pltargs)
    ax[1].semilogx(w,phase*180/np.pi,**pltargs)
bodeplot_set(ax,3,3)
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\freqplot.py:187: FutureWarning: 'Plot' keyword is deprecated in bode_plot; use 'plot'
  warnings.warn("'Plot' keyword is deprecated in bode_plot; use 'plot'",
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\freqplot.py:187: FutureWarning: 'Plot' keyword is deprecated in bode_plot; use 'plot'
  warnings.warn("'Plot' keyword is deprecated in bode_plot; use 'plot'",
C:\Users\h\AppData\Roaming\Python\Python39\site-packages\control\freqplot.py:187: FutureWarning: 'Plot' keyword is deprecated in bode_plot; use 'plot'
  warnings.warn("'Plot' keyword is deprecated in bode_plot; use 'plot'",

png

可以注意到,随着$\omega_n$的增大,其相移和衰减的程度都下降了,而在固有频率处,刚好相移为90°

1
2
P=tf([1,3],[1,3,2])
b=bode(P,logspace(-2,2))

​
png
​