- 这篇文章是现代控制理论课程设计的前半部分，主要研究了单自由度弹性关节系统的工程背景，动力学建模；控制系统方面包括状态方程建立，状态转移矩阵计算，零输入给定初值
MATLAB与Python Tkinter模型仿真，状态方程与传递函数之间的转换，系统能观性和能控性的讨论。 - 单自由度弹性关节示意图如下：
P.S. 这篇文章的原排版系统为LaTeX，在ShareLaTeX上可以导出PDF格式的报告，同时LaTeX Source也开放了公共阅读权限。
工程背景
工业机器人概述
工业机器人按照ISO8373定义\footnote{ISO Standard 8373:1994, Manipulating Industrial Robots – Vocabulary}，它是面向工业领域的多关节机械手或多自由度的机器人。工业机器人是自动执行工作的机器装置，是靠自身动力和控制能力来实现各种功能的一种机器。它可以接受人类指挥，也可以按照预先编排的程序运行，现代的工业机器人还可以根据人工智能技术制定的原则纲领行动。
工业机器人的重要特性是三维空间运动的空间机构，这也是其区别与数控机床的特征。机器人的机构本质上为空间多体系统。当考虑系统元件的弹性时，为多柔性体系统。机器人的主要机构有：关节、连杆、伺服电机、减速机等。
- 关节：机器人系统中将运动副称为关节。在关节之中，凡单独驱动的称为主动关节，反之称从动关节。工业机器人手臂的关节常为单自由度主动关节，即每一个关节均有一个驱动器驱动。
- 连杆：连杆机构是传递机械能的一种装置，通常是由刚体构件用转动副、移动副、球面副、球销副、圆柱副或螺旋副中的一种或几种联结而成的机械机构，因为上述联接副均属于低副，连杆机构也称为低副机构。通过不同的设计与计算，连杆机构可实现转动、直线移动、往复运动和平面或空间的复杂函数运动轨迹。
- 伺服电机：伺服电机是对用于使用伺服机构的电动机总称。伺服系统，就是依照指示命令动作所构成的控制装置，应用于马达的伺服控制，将感测器装在马达与控制对象机器上，侦测结果会返回伺服放大器与指令值做比较。伺服马达的动作特性是进行位置定位控制和动作速度控制，其主要特点为转速可以精确控制，速度控制范围广。
- 减速机：减速机一般用于低转速大扭矩的传动设备，把电动机，内燃机或其它高速运转的动力通过减速机的输入轴上的齿数少的齿轮啮合输出轴上的大齿轮来达到减速的目的，减速机是一种动力传达机构，利用齿轮的速度转换器，将马达的回转数减速到所要的回转数，并得到较大转矩的机构。
机器人柔性设计
机器人的柔顺性是实现未知约束环境下人机安全物理交互，开展复杂作业的 重要前提。针对机器人关节的顺应性控制问题，传统工业中普遍应用“刚性设计，控制实现安全”的模式，然而这种设计会造成操作者的不适，且有可能对操作者或机器人接触的物体造成损伤。为了解决这个问题，可以在驱动机构和负载端之间串联弹性元件，将负载输出和电机惯量隔离。在人机物理交互的场合中，可确保接触力控制在合理范围内，其实现模式为“设计实现安全，控制提高性能\footnote{浙江大学硕士学位论文：串联弹性关节控制与交互刚度辨识 吕铖杰 2015}”。但是弹性介质的引入增加了系统动态结构的复杂度，设计稳定和可靠的控制算法成为使用该类驱动器的挑战之一。
动力学建模
单自由度弹性关节的使用场景一般是低速大力矩情形下，往往电机和负载之间会安装减速器，为了分析方便，在理想状态下不考虑减速器摩擦与额外转动惯量的影响，以下为单自由度弹性关节示意图： 其中：\(J_m\)为电动机的转动惯量，\(J_l\)为负载的转动惯量，\(\beta_m\)为电动机的阻尼，可认为是转动摩擦与空气阻力等，\(\beta_l\)为负载的阻尼，\(\theta_m\)为电动机旋转的角度，\(\theta_l\)为负载旋转的角度，\(k\)为扭力弹簧常数，\(u\)为电动机的输入转矩。
首先对电动机进行受力分析：电动机输入转矩与弹簧扭转力矩，自身转动惯量产生的转矩，阻尼产生的转矩之和相平衡。
弹簧扭转力矩公式为：\(F_1=k\times\Delta\theta\)；
转动惯量产生的转矩为：\(F_2=J\times\omega\)；
阻尼产生的转矩为：\(F_3 = \beta\times\alpha\)；
综上，可导出电动机转动时的转矩平衡公式：
\[u = k(\theta_m-\theta_l)+J_m\ddot{\theta}_m+\beta_m \dot{\theta}_m\]
同理可导出负载转动时的转矩平衡公式：
\[k(\theta_l-\theta_m)=-J_l\ddot{\theta}_l-\beta_l\dot{\theta}_l\]
整理两式，单自由度弹性关节的模型为：
\[
\left\{ \begin{array}{ll}
J_l\ddot{\theta}_l+\beta_l\dot{\theta}_l+k(\theta_l-\theta_m)=0 \\
J_m\ddot{\theta}_m+\beta_m\dot{\theta}_m-k(\theta_l-\theta_m)=u
\end{array} \right.
\]
控制系统分析
微分方程转化状态方程
首先确定状态变量为：\(\theta_l,\dot{\theta}_l,\theta_m,\dot{\theta}_m\)；
根据前文导出的单自由度弹性关节模型，可导出系统的状态方程为：
\[
\left\{ \begin{array}{ll}
\dot{x_1} = x_2 \\
\dot{x_2} = -\frac{k}{J_l}x_1-\frac{\beta_l}{J_l}x_2+\frac{k}{J_l}x_3 \\
\dot{x_3} = x_4 \\
\dot{x_4} = \frac{k}{J_m}x_1-\frac{\beta_m}{J_m}x_4-\frac{k}{J_m}x_3+\frac{1}{J_m}u
\end{array} \right.
\] ## 状态转移矩阵 线性定常系统的状态空间表达式为：
\[
\left\{ \begin{array}{ll}
\dot{X} = AX+Bu \\
y = CX
\end{array} \right.
\]
其中：
\[A =
\begin{bmatrix}
0 & 1 & 0 & 0 \\
-\frac{k}{J_l} & -\frac{\beta_l}{J_l} & \frac{k}{J_l} & 0 \\
0 & 0 & 0 & 1 \\
\frac{k}{J_m} & 0 & -\frac{k}{J_m} & -\frac{\beta_m}{J_m}
\end{bmatrix};\qquad B =
\begin{bmatrix}
0 \\ 0 \\ 0 \\ \frac{1}{J_m}
\end{bmatrix}; \qquad C =
\begin{bmatrix}
1 & 0 & 0 & 0
\end{bmatrix}
\]
线性定常系统在没有控制作用时，由初始条件引起的运动为自由运动：对于线性定常系统，状态转移矩阵为： \[\Phi(t-t_0)=e^{A(t-t_0)}\]
MATLAB求解
代入系统已知参数，使用
|
|
或者
即可解出结果。
零输入给定初值状态方程仿真
单自由度弹性关节机器人系统仿真所需参数如下表所示:
| 物理量 | 数值 | 单位 |
|---|---|---|
| \(J_m\) | 1 | \(kg\cdot m^2\) |
| \(J_l\) | 1 | \(kg\cdot m^2\) |
| \(\beta_m\) | 0.1 | \(N\) |
| \(\beta_l\) | 0.1 | \(N\) |
| \(k\) | 60 | \(N/rad\) |
ODE45方法求解(基于MATLAB)
使用ode45求解高阶常微分方程需要做降阶处理，使用该系统的状态方程构造函数，代码如下：
|
|
其中，已知输入\(u\)为零。
代入矩阵参数，初值，时间求微分方程的数值解，并绘制\(\theta_l,\dot{\theta}_l,\theta_m,\dot{\theta}_m\)关于变量\(t\)的曲线，代码如下：
|
|
矩阵指数方法求解(基于MATLAB)
和上一种方案的矩阵参数，初值，时间相同，使用拉氏反变换法求解矩阵指数，并绘制\(\theta_l,\dot{\theta}_l,\theta_m,\dot{\theta}_m\)关于变量\(t\)的曲线。
矩阵指数的计算公式：
\[e^At=\mathcal{L}^{-1}[(sI-A)^{-1}]\]
MATLAB程序运行结果
运行以上MATLAB程序，可得到如下图所示的图像；当前给负载角度初值为\(0.1rad\),弹簧处于紧绷状态，系统零输入启动后，电机与负载产生相反方向的加速度，最终达到稳态。 
系统仿真动画(基于Python Tkinter)
首先使用Python scipy模块求解系统的状态方程，得到\(\theta_l,\dot{\theta}_l,\theta_m,\dot{\theta}_m\)关于变量\(t\)的变化情况。借由matplotlib模块绘制曲线。同时，使用Tkinter创建图形用户界面，并根据上述求解的值绘制动画，实现数据可视化。程序运行效果如下图所示：
代码功能：
单自由度弹性关节机器人系统的仿真程序，支持自定义输入，初值，系统参数；可以通过图形用户界面观察该系统的运行状态。
运行环境：
兼容Windows/macOS/Linux，需要Python环境：
- Python 2.7.11
- matplotlib \(\geq\) 1.5.1
- numpy \(\geq\) 1.12.1
- scipy \(\geq\) 0.19.0
代码说明:
首先导入需要使用的模块，相关代码块如下所示：
变量赋值部分，根据实际物理系统给转动惯量，阻尼，劲度系数等变量赋值，同时确定画布大小，确定代表电机/负载的矩形位置。
以下为数据初始化函数，规范数据格式：
以下为ode求解所需的降阶函数，输入为初始状态，时间，系统状态空间描述系数矩阵，输出重构为单维数组(Python odeinit解法不支持高阶输入)
以下为矩形旋转坐标运算函数，输入为矩形四点坐标，需要旋转的角度(弧度制)，旋转中心坐标；返回值为旋转后的矩形四点坐标。
以下为绘图函数，输入为几何元素canvas，ode解法得到的角度list，使用i索引list元素，无返回值。遍历解出的角度list即可画出旋转矩形；为实现动画效果，需使用定时任务，在该函数中使用了time模块的sleep()函数
以下为开始按钮的事件函数，输入为几何元素canvas，ode解法得到的角度list，Tk的界面对象。
以下为数据过程可视化函数，输入为ode解法得到的角度list。封装canvas元素，label元素(用以显示解释性文字)，buttom元素(开始绘图功能支持)。
以下为数据图表可视化函数，输入为ode解法得到的角度list。使用matplotlib模块中的plot函数绘制\(\theta_l,\dot{\theta}_l,\theta_m,\dot{\theta}_m\)关于变量\(t\)的曲线。
以下为主函数，初始化数据，求解模型的状态方程，并分别调用过程可视化函数，图表可视化函数。
状态方程转化为传递函数
传递函数可用来描述线性非时变系统的特性。一个连续线性非时变系统可通过以下方法将状态方程转化为传递函数。
首先对下式进行拉氏变换：
\[ \dot{X} = AX + Bu\] 可以得到：
\[sX(s) = AX(s) + Bu(s)\] 再针对\(X(s)\)进行化简可得到：
\[(sI-A)X(s)=Bu(s)\]
\[X(s) = (sI-A)^{-1}Bu(s)\]
用上式替换以下输出方程中的\(X(s)\)：
\[y(s) = CX(s)+Du(s)\]
结果如下：
\[y(s) = C((sI-A)^{-1}Bu(s))+Du(s)\]
由传递函数的定义\(G(s) = \frac{y(s)}{u(s)}\)：
可导出：
\[G(s) = C(sI-A)^{-1}B + D\]
传递函数 使用系统的状态方程的系数矩阵求解，求解出的结果为：
\[G(s) = \frac{k}{J_lJ_ms^4 + (B_lJ_m+ B_mJ_l)s^3 + (B_lB_m+J_lk+J_mk)s^2 + (B_lk+ B_mk)s }\]
代入已知系数，该系统的传递函数为：
\[G(s) = \frac{6000}{100s^4+20s^3+12001s^2+1200s}\]
能观性与能控性讨论
能控性讨论
定理： 系统\(\Sigma = (A,B)\)，即
\[\dot{X} = AX+Bu\]
\[y = CX+Du\]
状态完全能控的充分必要条件是其能控性矩阵
\[Q_k = \lbrack B \vdots AB \vdots \cdots \vdots A^{n-1}B \rbrack\]
满秩，即
\[rank\lbrack B \vdots AB \vdots \cdots \vdots A^{n-1}B \rbrack = n\]
由系统的状态方程系数矩阵\(A,B\)求出其能控性矩阵为：
\[Q_k =
\begin{bmatrix}
0 & 0 & 0 & \frac{k}{J_lJ_m} \\
0 & 0 & \frac{k}{J_lJ_m} & -\frac{B_lk}{J_mJ_l^2} - \frac{B_mk}{J_lJ_m^2} \\
0 & \frac{1}{J_m} & -\frac{B_m}{J_m^2} &-\frac{k}{J_m^2} - \frac{B_m^2}{J_m^3} \\
\frac{1}{J_m} & -\frac{B_m}{J_m^2} & -\frac{k}{J_m^2} -\frac{B_m^2}{J_m^3} & \frac{B_mk}{J_m^3} + \frac{B_mk}{J_m^3} - \frac{B_m^3}{J_m^4} \\
\end{bmatrix}
\]
该矩阵的行列式可以按下式求解：
\[det(Q_k) =
-\frac{k}{J_lJ_m} \times \frac{k}{J_lJ_m} \times
\begin{vmatrix}
0 & \frac{1}{J_m} \\ \frac{1}{J_m} & -\frac{B_m}{J_m^2} \\
\end{vmatrix}
=\frac{k^2}{J_l^2J_m^2}\frac{1}{Jm^2}
\]
由于\(J_l,J_m\neq 0\)，所以当且仅当\(k = 0\)时系统不能控。
能观性讨论
定理：系统\(\Sigma = (A,C)\)，即
\[\dot{X} = AX\]
\[y = CX\]
状态完全能观测的充分必要条件是其能能观测性矩阵
\[Q_g = \lbrack C \vdots A^TC^T \vdots \cdots \vdots (A^T)^{n-1}C^T \rbrack\]
满秩，即
\[Q_g =
\begin{bmatrix}
C \\
CA \\
\vdots \\
CA^{n-1} \\
\end{bmatrix}
\] 满秩。
由系统的状态方程系数矩阵\(A,C\)求出其能观测性矩阵为：
\[Q_g =
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
-\frac{k}{J_l} & -\frac{B_l}{J_l} & \frac{k}{J_l} & 0 \\
\frac{B_lk}{J_l^2} & \frac{B_l^2}{J_l^2} - \frac{k}{J_l} & -\frac{B_lk}{J_l^2} & \frac{k}{J_l}\\
\end{bmatrix}
\] 该矩阵的行列式可以按下式求解：
\[det(Q_g)=
1\times 1 \times
\begin{vmatrix}
\frac{k}{J_l} & 0 \\
-\frac{B_lk}{J_l^2} & \frac{k}{J_l} \\
\end{vmatrix}
= \frac{k^2}{J_l^2}
\] 由于\(J_l\neq 0\)，所以当且仅当\(k = 0\)时系统不能观测。
能控规范型和能观测规范型
单输入-单输出系统的能控规范型
单自由度弹性关节机器人系统的能控性矩阵： \[Q_k = \begin{bmatrix} 0 & 0 & 0 & \frac{k}{J_lJ_m} \\ 0 & 0 & \frac{k}{J_lJ_m} & -\frac{B_lk}{J_mJ_l^2} - \frac{B_mk}{J_lJ_m^2} \\ 0 & \frac{1}{J_m} & -\frac{B_m}{J_m^2} &-\frac{k}{J_m^2} - \frac{B_m^2}{J_m^3} \\ \frac{1}{J_m} & -\frac{B_m}{J_m^2} & -\frac{k}{J_m^2} -\frac{B_m^2}{J_m^3} & \frac{B_m k}{J_m^3} + \frac{B_m k}{J_m^3} - \frac{B_m^3}{J_m^4} \\ \end{bmatrix} \]为非奇异，则存在非奇异变换： \[\widehat{X}=PX\]
或
\[X=P^{-1}\widehat{X}\]
系统可化为能控规范型，即 \[P_1= \begin{bmatrix} 0 & 0 & 0 & 1 \\ \end{bmatrix}Q_k^{-1} \] 求得 \[P_1 = \begin{bmatrix} \frac{J_lJ_m}{k} & 0 & 0 & 0\\ \end{bmatrix} \] 变换矩阵为 \[ P = \begin{bmatrix} P_1 \\ P_1A \\ P_1A^2 \\ P_1A^3 \\ \end{bmatrix} = \begin{bmatrix} \frac{J_lJ_m}{k} & 0 & 0 & 0 \\ 0 & \frac{J_lJ_m}{k} & 0 & 0 \\ -J_m & -\frac{B_lJ_m}{k} & J_m & 0\\ \frac{B_l J_m}{J_l} & \frac{J_m B_L^2}{J_l k}-J_m & -\frac{B_l J_m}{J_l} & J_m \\ \end{bmatrix} \] 因此 \[ \widehat{A} = PAP^{-1}= \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & -\frac{k(B_l + B_m)}{J_l J_m} & -\frac{(J_l k + J_m k + B_l B_m)}{(J_l k + J_m k + B_l B_m)} & -\frac{B_l J_m + B_m J_l}{J_l*J_m} \\ \end{bmatrix} \qquad \widehat{B} = PB = \begin{bmatrix} 0 & 0 & 0 & 1 \\ \end{bmatrix}^T \] 故可将状态方程化为能控规范性 { \[\begin{array}{ll} \dot{\widehat{X}} = \widehat{A}\widehat{X}+\widehat{B}u \\ y = \widehat{C}\widehat{X} \end{array}\]
.
单输入-单输出系统的能观测规范型
单自由度弹性关节机器人系统的能观测性矩阵： \[Q_g^T =
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
-\frac{k}{J_l} & -\frac{B_l}{J_l} & \frac{k}{J_l} & 0 \\
\frac{B_lk}{J_l^2} & \frac{B_l^2}{J_l^2} - \frac{k}{J_l} & -\frac{B_lk}{J_l^2} & \frac{k}{J_l}\\
\end{bmatrix}
\]
为非奇异，则存在非奇异变换： \[X(t)=T\widehat{X}(t)\]
或 \[\widehat{X}(t)=T^{-1}X(t)\] 系统可化为能控规范型，即 \[T_1 = (Q_g^T)^{-1}
\begin{bmatrix}
0 & 0 & 0 & 1 \\
\end{bmatrix}^T
\]
变换矩阵为 \[T = \begin{bmatrix} T_1 & A T_1 & A^2 T_1 & A^3 T_1 \\ \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & -\frac{B_l}{J_l}-\frac{B_m}{J_m}\\ 0 & \frac{J_l}{k} & -\frac{B_m J_l}{J_m*k} & \frac{B_m^2 J_l}{J_m^2 k}-\frac{J_l}{J_m} \\ \frac{J_l}{k} & -\frac{B_m J_l}{J_m k} & \frac{B_m^2 J_l}{J_m^2 k}-\frac{J_l}{J_m} & \frac{k B_m J_l + k B_m J_m^2 - J_m B_m^3}{k J_m^3} \\ \end{bmatrix} \]
因此 \[\widehat{A} = T^{-1}AT = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & -\frac{k(B_l + B_m)}{J_l J_m} \\ 0 & 1 & 0 & -\frac{(J_l k + J_m k + B_l B_m}{J_l J_m} \\ 0 & 0 & 1 & -\frac{B_l J_m + B_m J_l}{J_l J_m}] \end{bmatrix} \qquad \widehat{C} = CT = \begin{bmatrix} 0 & 0 & 0 & 1 \\ \end{bmatrix} \] 故可将状态方程化为能观测规范性 { \[\begin{array}{ll} \dot{\widehat{X}} = \widehat{A}\widehat{X}+\widehat{B}u \\ y = \widehat{C}\widehat{X} \end{array}\]
.
控制系统设计
系统给出的参数为\(J_m=2,J_l=10,B_m=0.5,B_l=1,k=111\) ## 开环阶跃响应仿真 使用MATLAB将系统的状态空间方程转换为传递函数形式，并观察阶跃响应，代码如下：
Notes: 这里将输出矩阵\(C\)改为了四阶单位矩阵，这样可以观察系统各状态的阶跃响应情况。
阶跃响应曲线 
由阶跃响应曲线可以看出，电机和负载的转速基本保持同步，但二者相互影响；从角速度曲线上可以看出二者的角速度都存在锯齿，一段时间后保持稳定。电机和负载的转速均线性增加。
状态反馈控制器设计
根据性能指标要求确定系统主导极点
这是一个主导极点问题。在稳定的高阶系统中，对于其时间响应起主导作用的闭环极点为主导极点，远离原点的极点所对应的响应影响很小，离虚轴最近且邻近没有零点的极点，其对应的暂态响应分量衰减最慢，幅值也大。正是如此，将本例中讨论的四阶系统近似为二阶系统分析。系统要求的性能指标为超调量\(\sigma \le 5\%\);调节时间$ t_s 0.5s\(; 系统频宽\)BW 10$。
由超调量定义公式得到欠阻尼典型二阶系统的超调量为: \[\sigma \% = \frac{c(t_p)-1}{1}\times 100\% = e^{-\pi \xi/\sqrt{1-\xi^2}}\times 100\%\]
超调量\(\sigma \le 5\%\)，取超调量为\(1.52\%\),解得\(\xi = 0.8\)。
由调节时间公式: \[t_s \approx \frac{4}{\xi \omega_n}\]
调节时间间$ t_s 0.5s$，取调节时间为\(0.5s\),解得\(\omega_n = 10\)
由典型二阶系统的带宽频率公式： \[BW = \omega_n\sqrt{(1-2\xi^2)+\sqrt{(1-2\xi^2)^2+1}}\]
验证得出系统的带宽为\(8.7090\)。 由典型二阶欠阻尼系统的关系式： \[s_{1,2}=-\xi\omega_n=\pm j\omega_n\sqrt{1-\xi^2}\] \[s^2+2\xi\omega_n s+\omega_n^2=0\]
得出配置的闭环主导极点为： \[s_{1,2} = -8\pm 6j\]
远极点应选择使它和原点的距离大于\(5\mid s_1 \mid\)的点，所以确定希望的极点为：
{ \[\begin{array}{ll} s_1 = -8 + 6j \\ s_2 = -8 - 6j \\ s_3 = -50 \\ s_4 = -60 \\ \end{array}\]
.
确定状态反馈矩阵\(\widehat{K}\)
受控系统的特征多项式为 \[f(s)=s^4+0.35s^3+66.625s^2+8.325s;\qquad a_1=0.35\quad a_2=66.625\quad a_3=8.325\quad a_4 =0\] 由希望的极点构成的特征多项式为 \[f^*(s)=(s+50)(s+60)(s+8-6j)(s+8+6j) = s^4 + 126s^3 + 4860s^2 + 59000s + 300000\] \[a_1^* = 126 \quad a_2^* = 4860 \quad a_3^* = 59000 \quad a_4^* = 300000
\]
于是状态反馈矩阵\(\widehat{K}\)为 \[K=
\begin{bmatrix}
300000 & 58991.675 & 4793.375 & 125.65\\
\end{bmatrix}\] 若需将其转化为能控规范性反馈矩阵可按\(\widehat{K} = KP\)
确定输入放大系数\(L\)
已知闭环传递函数为 \[W_K(s)=\frac{5.55L}{s^4 + 126 s^3 + 4860 s^2 + 59000 s + 300000}\]
要求跟踪信号阶跃信号的误差\(e_p = 0\)，有
\[e_p = 0 = \lim\limits_{t\rightarrow{\infty}}[1-y(t)]=\lim\limits_{s\rightarrow{\infty}}s[\frac{1}{s}-\frac{W_K(s)}{s}]=\lim\limits_{s\rightarrow{\infty}}[1-W_K(s)] = \frac{300000-5.55L}{300000}=0\]
解得
\[L = 54055\]
MATLAB阶跃响应仿真
使用MATLAB进行极点配置及绘制单位阶跃响应曲线的代码如下：
绘制出的阶跃响应曲线如下图所示： 
并使用MATLAB计算超调量及调节时间，可用以下代码：
解出超调量为\(1.46\%\)，调节时间为\(0.37s\)，符合之前设置的性能指标要求。
带观测器状态反馈闭环系统
由于并非所有的状态变量都能够被测量到。因此，为了实现状态反馈控制律，就必须利用已知信息\((y,u)\)，通过模型来对系统的状态变量进行估计。 ### 观测器极点位置选择 观测器极点的选取通常使观测器的响应比系统的响应要快得多，一个经验法则是选择观测器的响应至少为系统响应的\(2\sim 5\)倍。观测器的极点位于所期望的闭环极点的左边，所以后者在响应中起主导作用。 根据以上规律，选择观测器极点为：{ \[\begin{array}{ll} s_1 = -70 \\ s_2 = -71 \\ s_3 = -72 \\ s_4 = -73 \\ \end{array}\]
.
确定观测器方程
设观测器增益为\(G = [g_2 g_1]^T\)，观测器期望的多项式为 \[a^*(s)=(s+70)(s+71)(s+72)(s+73)\] 观测器的特征多项式为： \[det[sI-(A-GC)]\] 可求观测器方程： \[\dot{\widehat{x}}=(A-G C)\widehat{x}+B u+G y\] 带观测器的状态反馈系统可以用分块矩阵表示： \[
\begin{pmat}[{.}]
\dot{x} \cr \-
\dot{x} - \dot{\widetilde{x}} \cr
\end{pmat}=
\begin{pmat}[{|}]
A-B K & BK \cr \-
0 & A-G C \cr
\end{pmat}
\begin{pmat}[{.}]
X \cr \-
X - \widetilde{X} \cr
\end{pmat}
+
\begin{pmat}[{.}]
B \cr \-
0 \cr
\end{pmat}
u
\] \[y=
\begin{pmat}[{|}]
C & 0 \cr
\end{pmat}
\begin{pmat}[{.}]
X \cr \-
X - \bar{X} \cr
\end{pmat}
\] ### MATLAB阶跃响应仿真 用MATLAB分别计算状态反馈矩阵和观测器反馈矩阵，构造分块矩阵并绘制阶跃响应曲线，代码如下：
带观测器状态反馈闭环系统阶跃响应曲线
这个响应结果几乎和之前没有采用观测器时的结果完全一致。这是因为我们假设的观测器模型与实际系统模型相同（包括有相同的初始条件）。所以，在最小的控制信号代价下，所有的设计要求都得到了满足，无需再次设计。
观测器期望极点比较性验证
在本例中，真实状态是无法观测的，但由于在MATLAB仿真中设置的观测器模型与系统模型完全相同，初值也相同，所以任意配置极点都能保证系统状态的准确估计。为了验证观测器极点位置对观测器逼近被估计状态速度的影响，可设置一个观测器初值误差，观察系统的响应情况。
由此。将观测器初值设置为\([1\quad0 \quad 0 \quad 0]\)，即负载转角初值未测准，误差\(1rad\)。 分别配置观测器极点为
.
使用以下代码配置初值，并分析阶跃响应情况：
得到的观测器响应曲线如下图所示： 
根据经典控制理论，我们知道闭环极点离开虚轴的距离越大，系统的调节时间越短。以上仿真的结果和理论是一致的，由图可以看出配置\(P_3\)极点，观测器的估计收敛速度最快。
系统的阶跃响应情况 
在初值误差相同的情况下，观测器闭环极点距离虚轴越远，响应速度越快，跟随性能越好。但是相对地，超调量也会增加。这同时也表征了观测器快速性和抗干扰性的矛盾。观测器的最大响应速度受到控制系统中的噪声和灵敏性的限制。为了保证系统有足够的抗扰性能，通常设计\(G\)使观测器比被估系统稍微快些就可以了。