0%

控制系统建模

1
2
3
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

描述控制系统

在这一部分,假设$u(t)$为系统某一时刻的输入,$y(t)$为系统某一时刻的输出$t$为时间

利用递推法即$y_n=\Sigma f(y_k)+g(u_n)$,我们可以得到一个齐次线性微分方程描述输出和输入之间的关系
$$\Sigma_k^n\frac{\mathbb{d}^k}{\mathbb{d}t^k}a_ky(t)=\Sigma_i^n\frac{\mathbb{d}^i}{\mathbb{d}t^i}b_ku(t)$$

对于一些特殊情况,例如系统的状态与其历史无关,那么系统可以被简化为$y(t)=b_0u(t)$
下面将对一个常见模型进行建模


垂直驱动机械臂模型

$$J\ddot{\theta}(t)+\mu\dot{\theta}(t)+Mgl\sin \theta(t)=u(t)$$
臂转动惯量为$J$,质量为$M$,重力加速度为$g$,考虑阻尼作用,整体外部扭矩被表示为$u$
但是这不是一个线性微分方程(由于$\sin$的存在),因此我们可以对$\theta$进行小角度假设,那么就可以得到一个线性微分方程

$$J\ddot{\theta}(t)+\mu\dot{\theta}(t)+Mgl\theta(t)=u(t)$$

控制工程所使用的模型描述

由于微分方程往往很高阶,因此直接分析会非常困难.
将微分方程转换为”传递函数模型”或者是”状态方程模型”

传递函数模型

传递函数模型形式$\mathcal{P}(s)=\frac{\Sigma^m_k b_ks^k}{\Sigma^n_i b_is^i}$
利用这种模型,可以将在时间域上的微积分操作转化为在$s$上的乘法操作,例如说$$\dot{y(t)}\to sy(s)$$
$\mathcal{P}(s)=\frac{y(s)}{u(s)}$其中$y(s)=\mathcal{L}[y(t)],u(s)=\mathcal{L}[u(t)]$

拉普拉斯变换$\mathcal{L}$,定义$$g(s)=\mathcal{L}[g(t)]:=\int_0^\infin g(t)e^{-st}\mathbb{d}t$$
对其进行$N$次积分相当于除以$s^N$,微分则相反

传递函数举例

对于之前的机械臂模型,我们可以将传递函数写作
$$\mathcal{P}(s)=\frac{1}{s^2J+\mu s+Mgl}$$
也可以利用Python中的(tf)函数来表示,例如我们现在有一个控制函数
$$\mathcal{P}(s)=\frac{1}{s^2+2s+3}$$
可以用以下代码表示

1
2
3
4
5
import control
np=[0,1]
dp=[1,2,3]
p=control.tf(np,dp)
print(p)
      1
-------------
s^2 + 2 s + 3

或者像这样$$\mathcal{P}(s)=\frac{s+2}{(s+1)(s+2)^2}$$

1
print(control.tf([1,2],[1,1])*control.tf([1],[1,2])**2)
        s + 2
---------------------
s^3 + 5 s^2 + 8 s + 4

可以利用(tfdata)重新获取分母和分子

1
2
[[nump]],[[denp]]=control.tfdata(p)
print(nump,denp)
[1] [1 2 3]

状态空间模型

状态空间模型可以将多元高阶微分方程表示为一阶微分方程的形式
通常可以编写为两个(两组)方程
$$
\dot{x}(t)=Ax(t)+Bu(t)\qquad (状态方程)
$$

$$
y(t)=Cx(t)+Du(t)\qquad (输出方程)
$$
需要注意的是,$ABCD$都可以是矩阵,$xyu$都可以是矢量,但是我们主要考虑单输入单输出情况(但是$x$可能还是矢量)

如何建模?
考虑$p=\frac{\mathbb{d}}{\mathbb{d}t}$,那么可以假设$p$的多项式$A(p)$和$B(p)$,于是我们就可以定义一个新变量$v$,有
$A(p)=u,y=B(p)v$
可以定义$x$向量为$v$的各阶导数,就可以写出原来的方程

机械臂的状态空间模型

下面以之前的机械臂为例,讨论一下状态空间模型的推导
由于机械臂运动状态由角度和角速度确定,因此$x=\left[\begin{matrix}\theta\\dot{\theta}\end{matrix}\right]$,然后其状态方程为
$$\dot{x}(t)=\left[\begin{matrix}\dot{\theta}\\ddot{\theta}\end{matrix}\right]=\left[\begin{matrix}0&1\-\frac{Mgl}{J}&-\frac{\mu}{J}\end{matrix}\right]x(t)+\left[\begin{matrix}0\\frac{1}{J}\end{matrix}\right]u(t)$$
输出方程为
$$y(t)=\left[\begin{matrix}1&0\end{matrix}\right]x(t)$$

状态空间模型的Python代码

可以使用函数ss(A,B,C,D)表示状态空间,非常方便,例如对于一组方程
$$\dot{x}(t)=\left[\begin{matrix}1&1&2\2&1&1\3&4&5\end{matrix}\right]x(t)+\left[\begin{matrix}2\0\1\end{matrix}\right]u(t)$$
$$y(t)=\left[\begin{matrix}1&1&0\end{matrix}\right]x(t)+u(t)$$

1
2
3
4
5
a="1 1 2;2 1 1;3 4 5"
b="2;0;1"
c="1 1 0"
d="1"
control.ss(a,b,c,d).A
array([[1., 1., 2.],
       [2., 1., 1.],
       [3., 4., 5.]])

框图

可以利用框图来描述控制系统

串联

将两个系统输入输出相连,可以得到一个串联的系统,例如将$y=\mathcal{S_1}u$与$z=\mathcal{S_2}y$串联后就得到$z=\mathcal{S_1}\cdot\mathcal{S_2}u$
系统可以表示为$$\mathcal{S=S_1\cdot S_2}$$
这一结果不满足交换律,框图如下所示
可以使用Python求得系统中串联的结果模型:

需要注意如果使用乘法,那么后作用的系统在前(参考矩阵),反之则是先作用的系统在前

并联

将输入分别通过两个过程再叠加(略)

反馈

当一个的输出再通过某些形式返回,就构成了反馈系统.例如在一个反馈系统中
$$y=\mathcal{S_1}r-\mathcal{S_2}y$$
那么这个反馈系统可以被化简为
$$y=\mathcal{S}r\qquad \mathcal{S}=\frac{\mathcal{S_1}}{1+\mathcal{S_1S_2}}$$

可以利用feedback实现反馈

1
2
3
4
s1=control.tf([0,1],[1,1])
s2=control.tf([1,1],[1,1,1])
s=control.series(s1,s2)
s

$$\frac{s + 1}{s^3 + 2 s^2 + 2 s + 1}$$