控制教程 —— 介绍篇：1.建模

介绍

从这篇文章开始，即进入正式的关于控制教程的介绍，当然介绍框架主要还是来自于密歇根大学的教程网站来，中间也会新增我自己的理解和补充。

正文

我们所要考察的系统是能够用外力作用来有目的地影响其运动行为的系统。这种系统叫做被控系统，或被控对象，简称对象。通常控制系统设计过程的第一步是要建立其数学模型，这些模型可以从物理定律或实验数据中得出。本节中，将用状态空间和传递函数来表示动态系统，并以一些机械和电气系统的例子来介绍建模的基本方法，并展示如何在MATLAB中生成这些模型来做进一步的分析。
本教程中使用的主要MATLAB命令是：ss， tf

动态系统

动态系统是根据固定规则随时间变化或演变的系统。对于许多物理系统，此规则可以表示为一组一阶微分方程：
$$\dot{\mathbf{x}}=\frac{d\mathbf{x}}{dt}=\mathbf{f}\left(\mathbf{x}(t)，\mathbf{u}(t)，t\right)$$在上面的等式中， $\mathbf{x}(t)$是状态变量，它是一组代表时间$t$的系统配置变量，例如，在简单的机械质量块-弹簧-阻尼器组成的系统中，两个状态变量可以是质量块的位置和速度。 $\mathbf{u}(t)$是在时刻$t$时由外部输入到系统的向量，而$\mathbf{f}$是一个(可能是非线性)函数，它表示在特定的时刻产生的状态向量对时间的导数(变化率)， $d\mathbf{x}/dt$。
当已知初始状态 $\mathbf{x}(t_0)$，以及时间 $t_0$ 和 $t_1$ 之间的所有输入 $\mathbf{u}(t)$，可通过积分上面的微分方程确定未来任何时刻的状态 $\mathbf{x}(t_1)$，尽管状态变量本身不是唯一的，可以有非常多个，但是为了能通过给定系统的“状态”，来预测系统未来的行为(求解状态方程)，需要确定一个最小数量的状态向量 $n$，它确定了状态空间的维度，通常会通过系统中独立储能元件的数量来确定 $n$。
前面的微分方程给出的状态关系比较一般化，可用于描述各种不同的系统，不幸的是，这可能导致分析起来很困难。通常采用两种简化方法来便于处理。

• 首先，如果函数 $\mathbf{f}$与时间无关，也就是 $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x}，\mathbf{u})$，这样的系统被称为时不变系统，这也是一个比较合理的假设，因为基本的物理定律本身并不依赖于时间。对于时不变系统，函数 $\mathbf{f}$的参数或系数是恒定的，但是状态变量 $\mathbf{x}(t)$和输入控制量 $\mathbf{u}(t)$可能仍然与时间有关。
• 第二个常见的是假设系统是线性的，实际上，几乎每个物理系统都是非线性的。换句话说，$\mathbf{f}$通常是状态和输入的某些复杂函数，这些非线性以多种不同的方式出现，控制系统中最常见的一种是“饱和”，指系统的某个量达到了物理极限。幸运的是，在足够小的工作范围内(理想曲线附近的切线)，大多数系统的动力学近似线性，在这种情况下，一阶微分方程组可以表示为矩阵方程，即 $\dot{\mathbf{x}}=A\mathbf{x}+B\mathbf{u}$。

在数字计算机问世之前(以及随后的大部分时间)，分析线性时不变(LTI)系统是切实可行的，因此，控制理论的大多数结果都是基于这些假设，幸运的是，正如我们看到的那样，这些结果被证明是非常有效的，并且使用LTI技术已经解决了许多重大的工程挑战。实际上，反馈控制系统的真正厉害之处在于它们在不可避免的不确定性存在的情况下建模仍能正常工作(鲁棒)。

状态空间表示

对于连续线性时不变(LTI)系统，标准状态空间表示如下：
$$\dot{\mathbf{x}}=A\mathbf{x}+B\mathbf{u}$$$$\mathbf{y}=C\mathbf{x}+D\mathbf{u}$$其中 $\mathbf{x}$ 是状态向量($n\times 1$)，$\dot{\mathbf{x}}$是状态向量对时间的导数($n\times 1$)，$\mathbf{u}$是输入或控制向量($p\times 1$)，$\mathbf{y}$是输出向量($q\times 1$)，$A$是系统矩阵($n\times n$)，$B$是输入矩阵($n\times p$)，$C$是输出矩阵($q\times n$)，$D$是前馈矩阵($q\times p$)。
输入矩阵有时非必要，因为通常存在有些状态变量无法直接观察或者不感兴趣，输出矩阵 $C$用于指定哪些状态变量(或组合)可供控制器使用。同样，通常情况下，输出并不直接取决于输入(仅通过状态变量)，在这种情况下，$D$是零矩阵。
状态空间表示(也称为时域表示)可以轻松处理**多输入/多输出(MIMO)**系统，初始条件非零的系统和非线性系统。因此，状态空间表示法在现代控制理论中广泛使用。

传递函数表示

LTI系统具有极其重要的特性，即如果系统的输入为正弦波，则输出也将为正弦波，其频率与输入频率相同，但幅值和相位可能不同。这些幅度和相位差是频率的函数，这就是所谓的系统频率响应。
使用拉普拉斯变换，可以将系统的时域表示形式转换为频域输入/输出表示形式，即传递函数。这样，它还可以将控制系统的微分方程转换成更易于分析的代数方程。
时域函数$f(t)$的拉普拉斯变换定义如下：
F ( s ) = L { f ( t ) } = ∫ 0 ∞ e − s t f ( t ) d t F(s) = \mathcal{L}\{f(t)\} = \int_0^\infty e^{-st}f(t)dt F(s)=L{f(t)}=∫0∞​e−stf(t)dt其中参数$s=\sigma +j\omega$是复数频率变量，在实际中，很少需要直接计算拉普拉斯变换，而是通过查表来转换，可以查看Laplace转换表(wiki百科)
$n$阶导数的拉普拉斯变换特别重要：
$$\mathcal{L}\left\{\frac{d^{n}f}{d{t}^n}\right\}=s^{n}F(s)-s^{n-1}f(0)-s^{n-2}\dot{f}(0)-...-f^{(n-1)}(0)$$频域方法最常用于分析LTI单输入/单输出(SISO)系统，例如微分方程的系数为常数的那些，如下所示：
$$a_n\frac{d^{n}y}{d{t}^n}+...+a_1\frac{dy}{dt}+a_{0}y(t)=b_m\frac{d^{m}u}{d{t}^m}+...+b_1\frac{du}{dt}+b_{0}u(t)$$该方程式的拉普拉斯变换如下：
$$a_n{s}^{n}Y(s)+...+a_{1}sY(s)+a_{0}Y(s)=b_m{s}^{m}U(s)+...+b_{1}sU(s)+b_{0}U(s)$$其中$Y(s)$ 和 $U(s)$分别是$y(t)$和$u(t)$的拉普拉斯变换，请注意，在查找传递函数时，我们始终假定每个初始条件$y(0)$，$\dot{y}(0)$，$u(0)$均为0。因此$$G(s)=\frac{Y(s)}{U(s)}=\frac{b_m{s}^m+b_{m-1}{s}^{m-1}+...+b_{1}s+b_0}{a_n{s}^n+a_{n-1}{s}^{n-1}+...+a_{1}s+a_0}$$将传递函数的分子和分母分解为零极点增益形式是很有用的：
$$G(s)=\frac{N(s)}{D(s)}=K\frac{(s-z_1)(s-z_2)...(s-z_{m-1})(s-z_m)}{(s-p_1)(s-p_2)...(s-p_{n-1})(s-p_n)}$$传递函数的零点$z_1,\dots ,z_m$，是分子多项式的根，也就是使 $N(s)=0$的$s$。传递函数的极点$p_1,\dots ,p_n$是分母多项式的根。也就是使$D(s)=0$的$s$，零点和极点都有可能是复数值(具有实部和虚部)，系统增益为 $K=b_m/a_n$。
请注意，我们还可以直接根据状态空间表示来确定传递函数，如下：
$$G(s)=\frac{Y(s)}{U(s)}=C(s\mathbf{I}-A{)}^{-1}B+D$$

机械系统

牛顿运动定律构成了分析机械系统的基础，牛顿的第二定律，可以表示为：
$$\Sigma \mathbf{F}=m\mathbf{a}=m\frac{d^2\mathbf{x}}{d{t}^2}$$指出作用在物体上的合力和等于其质量与加速度的乘积。为了我们的目的，牛顿第三定律指出，如果两个物体接触，那么它们将承受相同大小的接触力，只是作用方向相反。
在使用该方程时，最好构造一个显示所有作用力的受力图(free-body diagram,FBD)
例：
质量块弹簧阻尼系统

该系统的受力图如下：

弹簧的力与质量块的位移 $x$ 成正比，粘性阻力与质量块的速度 $v=\dot{x}$ 成正比，两种力都与质量块的运动相反，因此显示为 $x$ 的负方向。还应注意 $x=0$ 对应于弹簧未拉伸时的质量块位置。
现在，我们对力进行求和并在每个方向上应用牛顿第二定律，在这种情况下，没有作用在 $y$ 方向上的力，但是在 $x$ 方向上，我们有：
$$\Sigma F_x=F(t)-b\dot{x}-kx=m\ddot{x}$$该方程式称为控制方程式，完全表征了系统的动态状态。稍后，我们将看到如何使用它来计算系统对任何外部输入$F(t)$的响应，以及分析系统属性(例如稳定性和性能)。
为了确定质量块-弹簧-阻尼系统的状态空间表示，我们必须将二阶控制方程式简化为一组两个一阶微分方程式。为此，我们选择位置和速度作为状态变量。
$$\mathbf{x}=\left[\begin{array}{c}x\\ \dot{x}\end{array}\right]$$位置变量表征存储在弹簧中的势能，而速度变量表征存储在质量块中的动能。阻尼器即消耗能量，也不存储能量。通常，在选择状态变量时，考虑哪些变量会表征系统中存储的能量会很有帮助。
在这种情况下，状态方程：
$$\dot{\mathbf{x}}=\left[\begin{array}{c}\dot{x}\\ \ddot{x}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -\frac{k}{m}& -\frac{b}{m}\end{array}\right]\left[\begin{array}{c}x\\ \dot{x}\end{array}\right]+\left[\begin{array}{c}0\\ \frac{1}{m}\end{array}\right]F(t)$$如果我们对质量块的位置感兴趣，则输出方程可以表示为：
$$y=\left[\begin{array}{cc}1& 0\end{array}\right]\left[\begin{array}{c}x\\ \dot{x}\end{array}\right]$$
在MATLAB中输入状态空间模型
现在，我们将演示如何将上面导出的方程式输入到MATLAB的m文件中，让我们为每个变量分配以下数值

变量	符号	大小
m	质量	1.0 kg
k	弹簧常数	1.0 N/m
b	阻尼常数	0.2 Ns/m
F	输入力	1.0 N

创建一个新的m文件并输入一下命令

m = 1;
k = 1;
b = 0.2;
F = 1;

A = [0 1; -k/m -b/m];
B = [0 1/m]';
C = [1 0];
D = [0];

sys = ss(A,B,C,D)

假设初始条件为零，则该系统的拉普拉斯变换为：
$$m{s}^{2}X(s)+bsX(s)+kX(s)=F(s)$$因此，从输入力到位移输出的传递函数为：
$$\frac{X(s)}{F(s)}=\frac{1}{m{s}^2+bs+k}$$在MATLAB中输入传递函数模型
现在，我们将演示如何在MATLAB中创建以上导出的传递函数模型。在定义了系统参数的m文件中输入以下命令：

s = tf('s');
sys = 1/(m*s^2+b*s+k)

请注意，我们在这里使用了符号s来定义传递函数模型，我们建议大多数时候使用这种方法，但是在某些情况下，例如在较旧的版本的MATLAB中或与SIMULINK接口时，可能需要直接使用分子和分母多项式系数来定义传递函数模型，在这种情况下，请使用以下命令：

num = [1];
den = [m b k];
sys = tf(num,den)

电气系统

就像牛顿的机械系统定律一样，基尔霍夫的电路定律是对电气系统进行建模的基本分析工具，基尔霍夫(Kirchoff)的电流定律(KCL)规定，进入电路中某个节点的电流之和必须等于离开该节点的电流之和。基尔霍夫(Kirchoff)的电压定律(KVL)规定，电路中任何闭环周围的电压差之和为零。当施加KVL时，通常将源级电压取为正，将负载电压取为负。
例：
RLC电路
现在，我们将考虑三个无源电子元件的简单串联组合：电阻器，电感器和电容器，称为RLC电路。

由于该电路是一个单回路，因此每个节点只有一个输入和输出，因此，KCL的应用只是表明在任何给定时间 $i(t)$，整个电路的电流是相同的。现在在回路周围应用KVL并使用图中所示的符号约定，我们得出以下控制方程式：
$$V(t)-Ri-L\frac{di}{dt}-\frac{1}{C}\int idt=0$$我们注意到，RLC电路的控制方程与质量块-弹簧-阻尼器机械系统具有相似的形式。特别地，它们都是二阶系统，其中电荷(电流的积分)对应于位移，电感对应于质量，电阻对应与粘性阻尼，电容对应于弹簧刚度。在概念上，这些类比在理解动力学系统的行为方面非常有用。
通过选择电容器上的电荷和通过电路(电感器)的电流作为状态变量，可以找到状态空间表示形式：
$$\mathbf{x}=\left[\begin{array}{c}q\\ i\end{array}\right]$$其中$$q=\int idt$$因此，状态方程为：
$$\dot{\mathbf{x}}=\left[\begin{array}{c}i\\ \frac{di}{dt}\end{array}\right]=\left[\begin{array}{cc}0& 1\\ -\frac{1}{LC}& -\frac{R}{L}\end{array}\right]\left[\begin{array}{c}q\\ i\end{array}\right]+\left[\begin{array}{c}0\\ \frac{1}{L}\end{array}\right]V(t)$$我们选择电流作为输出，如下所示：
$$y=\left[\begin{array}{cc}0& 1\end{array}\right]\left[\begin{array}{c}q\\ i\end{array}\right]$$传递函数表示可以通过采用拉普拉斯变换来获得，就像我们对质量块弹簧阻尼器所做的那样，或者从状态空间方程获得：
$$\frac{I(s)}{V(s)}=C(s\mathbf{I}-A{)}^{-1}B+D=\left[\begin{array}{cc}0& 1\end{array}\right]{\left(s\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]-\left[\begin{array}{cc}0& 1\\ -\frac{1}{LC}& -\frac{R}{L}\end{array}\right]\right)}^{-1}\left[\begin{array}{c}0\\ \frac{1}{L}\end{array}\right]$$$$\Rightarrow \frac{I(s)}{V(s)}=\frac{s}{L{s}^2+Rs+\frac{1}{C}}$$可以使用与上述质量块弹簧阻尼器系统相同的步骤，将RLC状态空间和传递函数模型输入MATLAB。

系统辨识

在本节中，我们已经看到了如何使用基本物理原理对系统进行建模，但是，通常这是不可能的，因为系统的参数不确定，或者根本不了解其基本过程，在这些情况下，我们必须依靠实验测量和统计技术来开发系统模型，这一过程称为系统辨识。
可以使用时域或频域数据执行系统辨识，有关更多更多信息，可参考http://ctms.engin.umich.edu/CTMS/index.php?aux=Extras_Identification，同时今后会详细介绍

系统转换

MATLAB中的大多数运算可以在传递函数，状态空间模型或零极点增益形式上执行，此外，如果需要其他表示形式，则可以在这些形式之间进行转换。可参考http://ctms.engin.umich.edu/CTMS/index.php?aux=Extras_Conversions