转载自–如何入门现代控制理论 - J Pan的文章 - 知乎
https://zhuanlan.zhihu.com/p/57051153
一、什么是特征值和特征向量
二、什么是特征值和特征函数
三、什么是状态方程
四、如何求解状态方程
五、可控性和可观性
六、如何用状态反馈进行极点配置
七、如何进行跟踪控制
八、全状态观测器
九、缩减状态观测器
十、小结

一、什么是特征值和特征向量

相信大家都学过矩阵，这是大学最基本的课程，矩阵里面有一个核心的概念，叫特征值和特征向量，怎么定义的呢？大概就长这个样子：

Av=λvAv=\lambda vAv=λv

其中 AAA是一个 n×nn\times nn×n 的矩阵， vvv 是一个 n×1n\times 1n×1 的向量， λ\lambdaλ 是一个标量。可不要小看这个方程，它蕴藏着巨大的宝藏，要深入理解这式子并不简单。我们先从方程的左半部分 AvAvAv 说起，即什么是矩阵的乘法，举个简单的例子，假设我们现在有一个向量在 xxx 轴上，
v=[10]v=\begin{bmatrix} 1\\ 0 \end{bmatrix}v=[10]

同时有一个矩阵

A=[0110]A=\begin{bmatrix} 0 &1\\ 1 &0 \end{bmatrix}A=[0110]

那两者乘积是什么呢？很容易得到;

Av=[0110][10]=[01]=bAv=\begin{bmatrix} 0 &1\\ 1 &0 \end{bmatrix} \begin{bmatrix} 1\\ 0 \end{bmatrix} =\begin{bmatrix} 0\\ 1 \end{bmatrix}=bAv=[0110][10]=[01]=b

可以看到，在 xxx 轴上的向量 v=[10]Tv=[1 \quad 0]^Tv=[10]T通过矩阵 AAA 的乘积，变成了在 yyy 轴上的向量 b=[01]Tb=[0 \quad 1]^Tb=[01]T 。

因此，我们可以得出第一个结论：矩阵代表着某种变换，可以将空间中的一个向量变成另外一个向量。

好，现在我们要从另外一个角度来看，我们知道两个向量 a=[a1a2]a=[a_1 \quad a_2 ]a=[a1a2] ，b=[b1b2]b=[b_1 \quad b_2 ]b=[b1b2] 的内积定义为：

a⋅b=a1b1+a2b2a\cdot b=a_1b_1+a_2b_2a⋅b=a1b1+a2b2

矩阵的乘法是怎么定义的呢？

[a11a12a21a22][b1b2]=[a11b1+a12b2a21b1+a22b2]\begin{bmatrix} a_{11} &a_{12}\\ a_{21} &a_{22} \end{bmatrix} \begin{bmatrix} b_1\\ b_2 \end{bmatrix} =\begin{bmatrix} a_{11}b_1+a_{12}b_2\\ a_{21}b_1+a_{22}b_2 \end{bmatrix}[a11a21a12a22][b1b2]=[a11b1+a12b2a21b1+a22b2]

可以看到，矩阵的乘法其实就是内积，而内积代表什么？——投影啊！所以，我们就可以得到第二个结论：矩阵可以看成是一个坐标系，矩阵和向量的乘积，就是该向量在坐标系的投影。

矩阵一会是坐标系，一会又变成了变换，那到底是怎么一回事？——这就有点像波粒二象性了，光到底是粒子还是波？——物理学家是怎么解决这个问题的？——在经典物理学里面，光既是粒子，又是波，看你怎么观测它。矩阵也一样，站在不同的角度，它的作用就不一样。

怎么说？——假如我们现在有一个向量 vvv ，我们想变成向量 bbb ，怎么办？——两种方式，第一：坐标系不变，把向量 vvv 转换到 bbb ，对应第一条矩阵的解释；第二：向量不变，把坐标系变换了，向量vvv 在另一个坐标系下就变成了 bbb 。两个完全等价的，因为矩阵描述的变换，而变换是相对的。

好了，我们费了好大的劲，终于把 Av=λvAv=\lambda vAv=λv 的左半部分说完了，现在我们再看一下完整的方程又代表什么意思，有了前面的铺垫，我们解释起来就简单多了。这个方程啥意思呢？——向量 vvv 通过 AAA 的变换，幅值变化了 λ\lambdaλ 倍，方向未发生变化；换个说法，存在向量 vvv ，其在坐标系 AAA 中的投影，方向和 vvv 相同，幅值变化了 λ\lambdaλ 倍。也就是说，我们可以用很多 vˉ\bar vvˉ 来对 AAA 进行测试，结果发现，坐标系 AAA 不再是一团乱麻，在某些方向上出现了一些特殊的结果。比如在向量 vvv 的方向，出现了方向和数值的解耦，这就给我们提供了某种坐标系 AAA 的信息，这些信息是 AAA 的内部固有结构决定的，所以称之为特征值和特征向量。

那我们该怎样利用这些信息呢？可以证明，特征值和特征向量都不止一个，不妨假设它的特征向量别为 v=[v1,v2,⋯ ,vn]\bm{v}=[v_1, v_2 , \cdots ,v_n]v=[v1,v2,⋯,vn] ，

对应的特征值分别为：

λ=[λ10⋯00λ2⋯0⋮⋮⋮⋮00⋯λn]\bm{\lambda}=\begin{bmatrix} \lambda_1 &0& \cdots &0\\ 0 &\lambda_2 & \cdots & 0\\ \vdots& \vdots & \vdots &\vdots\\ 0 & 0 & \cdots & \lambda_n \end{bmatrix}λ=⎣⎢⎢⎢⎡λ10⋮00λ2⋮0⋯⋯⋮⋯00⋮λn⎦⎥⎥⎥⎤

则有：

Av=vλA \bm{v}=\bm{v}\bm{\lambda}Av=vλ

两边左边同时乘以 v−1\bm{v}^{-1}v−1 ，就可以得到：

v−1Av=λ\bm{v}^{-1}{}A \bm{v}=\bm{\lambda}v−1Av=λ

也就是说，通过特征向量组成的矩阵，我们顺藤摸瓜，顺便就把矩阵 AAA 变成了一个对角矩阵，而对角矩阵是一个性质非常好的矩阵，尤其适合当坐标系使用，比如最常用的三维欧式空间的坐标系就可以表示为：

I=[100010001]I=\begin{bmatrix} 1 &0 &0\\ 0 &1 & 0\\ 0 & 0 & 1 \end{bmatrix}I=⎣⎡100010001⎦⎤

总结一下：矩阵为我们提供了一种描述“变化”的方法，而特征变换则是利用矩阵本身具有的固有特性，将矩阵代表的坐标系转换到更简单的坐标系(或者是基)来进行研究，从而使“变化”这种现象变得更加简单。

二、什么是特征值和特征函数

前面我们说了，矩阵的特征向量有如下特性：

Av=λvAv=\lambda vAv=λv

矩阵 AAA 代表某种转换，把变量 vvv 变成 λv\lambda vλv 。现在我们从一个更数学的角度来看一下。假设AAA代表某种映射， vvv 也不再是向量，而是函数，即

A(x)=λxA(x)=\lambda xA(x)=λx

映射可以有很多种，不同的映射规则可以得到不同的结果，比如说，映射 AAA 代表微分算子 ddt\frac{d}{dt}dtd ，那会是怎么样呢？我们不妨简单推导一下看看：

dx(t)dt=λx(t)\frac{dx(t)}{dt}=\lambda x(t)dtdx(t)=λx(t)

两边同时乘以 dtx(t)\frac{dt}{x(t)}x(t)dt 可以得到：

dx(t)x(t)=λdt\frac{dx(t)}{x(t)}=\lambda dtx(t)dx(t)=λdt

两边积分就可以得到：

x(t)=x(0)eλtx(t)=x(0)e^{\lambda t}x(t)=x(0)eλt

也就是说，对于微分算子 ddt\frac{d}{dt}dtd ， eλte^{\lambda t}eλt 是其天然的特征函数，因为对于指数函数来说， ddt(eλt)=λeλt\frac{d}{dt}{\big(e^{\lambda t}}\big)=\lambda e^{\lambda t}dtd(eλt)=λeλt ，对于微分运算来说，运算的结果只是幅值变化了。当然，我还可以对 x(t)x(t)x(t) 的数量进行扩展，比如：

x(t)=[x1(t)x2(t)⋮xn(t)]\bm{x(t)}=\begin{bmatrix} x_1(t)\\ x_2(t)\\ \vdots \\ x_n(t) \end{bmatrix}x(t)=⎣⎢⎢⎢⎡x1(t)x2(t)⋮xn(t)⎦⎥⎥⎥⎤

则得到：

x˙(t)=λx(t)\bm{\dot x}(t)=\bm{\lambda} \bm{x}(t)x˙(t)=λx(t)

写成展开的形式更容易看一些：

ddt[x1(t)x2(t)⋮xn(t)]=[λ10⋯00λ2⋯0⋮⋮⋮⋮00⋯λn][x1(t)x2(t)⋮xn(t)]\frac{d}{dt}\begin{bmatrix} x_1(t)\\ x_2(t)\\ \vdots \\ x_n(t) \end{bmatrix} =\begin{bmatrix} \lambda_1 &0& \cdots &0\\ 0 &\lambda_2 & \cdots & 0\\ \vdots& \vdots & \vdots &\vdots\\ 0 & 0 & \cdots & \lambda_n \end{bmatrix}\begin{bmatrix} x_1(t)\\ x_2(t)\\ \vdots \\ x_n(t) \end{bmatrix}dtd⎣⎢⎢⎢⎡x1(t)x2(t)⋮xn(t)⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡λ10⋮00λ2⋮0⋯⋯⋮⋯00⋮λn⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡x1(t)x2(t)⋮xn(t)⎦⎥⎥⎥⎤

很容易得到，微分方程组的解为：

[x1(t)x2(t)⋮xn(t)]=[x1(0)eλ1tx2(0)eλ2t⋮xn(0)eλnt]\begin{bmatrix} x_1(t)\\ x_2(t)\\ \vdots \\ x_n(t) \end{bmatrix} =\begin{bmatrix} x_1(0)e^{\lambda_1t}\\ x_2(0)e^{\lambda_2t}\\ \vdots \\ x_n(0)e^{\lambda_nt} \end{bmatrix}⎣⎢⎢⎢⎡x1(t)x2(t)⋮xn(t)⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡x1(0)eλ1tx2(0)eλ2t⋮xn(0)eλnt⎦⎥⎥⎥⎤

也就说，假设所有满足 x˙(t)=λx(t)\bm{\dot x}(t)=\bm{\lambda} \bm{x}(t)x˙(t)=λx(t) 的函数组成一个空间，则这个空间的变量 x(t)\bm{x}(t)x(t) 是由仅由特征函数 [eλ1teλ2t⋯eλnt]\begin{bmatrix} e^{\lambda_1t}& e^{\lambda_2t}& \cdots &e^{\lambda_nt} \end{bmatrix}[eλ1teλ2t⋯eλnt] 叠加组成。换句话说， [eλ1teλ2t⋯eλnt]\begin{bmatrix} e^{\lambda_1t}& e^{\lambda_2t}& \cdots &e^{\lambda_nt} \end{bmatrix}[eλ1teλ2t⋯eλnt] 是 x˙(t)=λx(t)\bm{\dot x}(t)=\bm{\lambda} \bm{x}(t)x˙(t)=λx(t) 解空间的“基”或者“坐标系”，而[x1(0)x2(0)⋯xn(0)]\begin{bmatrix} {x_1(0)}& {x_2(0)}& \cdots &{x_n(0)} \end{bmatrix}[x1(0)x2(0)⋯xn(0)]则可以看成是坐标值。

特征函数和特征向量一样，就是利用“变换”本身的特性，寻找一种更简单的“坐标（或者说基）”来解决问题。而对于微分这种运算来说，指数函数 eλte^{\lambda t}eλt 就是最简单的一种“坐标”。

当然，对于坐标而言， x(t)\bm{x(t)}x(t) 只要求独立，不要求一定正交，因此矩阵 λ\bm{\lambda}λ 不再是对角矩阵，可以是更一般的矩阵，我们不妨用 AAA 表示，这样我们就得到：

x˙(t)=Ax(t)\bm{\dot x}(t)=A \bm{x}(t)x˙(t)=Ax(t)

眼尖的童鞋可能已经看到，这不就是状态方程吗？——没错！

三、什么是状态方程

前面我们已经基本得到了状态方程了，可是仔细一看呢，这个方程和输入貌似没有关系，那可以看成是零输入，现在我们要把输入加上：

x˙(t)=Ax(t)+Bu(t)\bm{\dot x}(t)=A \bm{x}(t)+B\bm{u(t)}x˙(t)=Ax(t)+Bu(t)

我们说了， [x1(t),x2(t),⋯ ,xn(t)]′[x_1(t),x_2(t),\cdots,x_n(t)]'[x1(t),x2(t),⋯,xn(t)]′ 可以看成是变量，导数 x˙(t)\bm{\dot x}(t)x˙(t)代表变化，整个式子放在一起，就是系统在在输入 u(t)\bm{u}(t)u(t) 的作用下，变量 x(t)\bm{x}(t)x(t) 在坐标系 [eλ1teλ2t⋯eλnt]\begin{bmatrix} e^{\lambda_1t}& e^{\lambda_2t}& \cdots &e^{\lambda_nt} \end{bmatrix}[eλ1teλ2t⋯eλnt] 中是如何变化的。

当然我们还关心在坐标系下的输出是啥，可以这么定义：

y(t)=Cx(t)+Du(t)\bm{y}(t)=C\bm{x}(t)+D\bm{u}(t)y(t)=Cx(t)+Du(t)

称之为输出方程，矩阵 CCC 可以看成是坐标值， Du(t)D\bm{u}(t)Du(t) 可以看成是前馈，即输入直接到输出了。画成框图的形式如下：

很多系统没有前馈，也就是 D=0D=0D=0 ，这个时候可简化为：

这种通过状态方程来描述系统的方式，我们一般称现代控制理论，那和经典控制理论有什么区别呢？首先先说相同点：那就是这两者应该是等价的，因为都是对自然现象的描述。区别就是现代控制理论更精细。啥意思呢？——打个比方，我们现在面前有一栋楼，经典控制理论会说：这栋楼有9层，在1楼有一个入口，一个出口，入口流量是xx，出口流量是xx，即SISO（single input single output）；而现代控制理论会怎么说呢？它会说，这栋楼有9层，每层都有一个入口和出口，第一层的出、入口流量是xx，第二层出、入口的流量是xx，…第9层出、入口的流量是xx，即MIMO（multiple input multiple output）。也就是说，现代控制理论的状态方程，用矩阵刻画了系统的内部结构，将系统的变化轨迹投影到更简单的坐标系上，从而可以在各个维度上进行描述。

先以简单的直流电机为例，一个典型的直流电机的数学模型如下：

其中 KmK_mKm 为力矩系数， KbK_bKb 为反电势系数， LLL 为线圈电感， RRR 为线圈电阻， JJJ 为转子转动惯量， KfK_fKf 为转子阻尼系数， TdT_dTd 为负载， ω\omegaω 为转子速度， θ\thetaθ 为转子位置， VaV_aVa 为端部电压，这就是理想直流电机的数学模型，注意这是开环的电机模型。

电压平衡方程式为： Va(t)=Ri(t)+Ldi(t)dt+Vemf(t)V_a(t)=Ri(t)+L\frac{di(t)}{dt}+V_{emf}(t)Va(t)=Ri(t)+Ldtdi(t)+Vemf(t)

反电势方程为： Vemf(t)=Kbθ˙(t)=Kbω(t)V_{emf}(t)=K_b\dot{\theta}(t)=K_b\omega (t)Vemf(t)=Kbθ˙(t)=Kbω(t)

力矩方程为： Ta(t)=Km⋅i(t)T_a(t)=K_m \cdot i(t)Ta(t)=Km⋅i(t)

转子的力矩平衡方程为： Jθ¨(t)=Ta(t)−Td(t)−Kfθ˙(t)J\ddot{\theta}(t)=T_a(t)-T_d(t)-K_f\dot{\theta}(t)Jθ¨(t)=Ta(t)−Td(t)−Kfθ˙(t)

如果都采用SI单位制，则 Km=Kb=KK_m=K_b=KKm=Kb=K ，假定状态变量为转角 θ(t)\theta(t)θ(t) （代表位置），转速 ω(t)\omega(t)ω(t) （代表速度）以及电流 i(t)i(t)i(t) （代表加速度），则于是可以获得直流电机的状态方程为：

[θ˙(t)ω˙(t)i˙(t)]⎵x˙(t)=[0100−KfJKJ0−KL−RL]⎵A[θ(t)ω(t)i(t)]⎵x(t)+[001L]⎵Bu(t)\underbrace{\begin{bmatrix} \dot\theta(t) \\ \dot\omega(t) \\ \dot i(t)\end{bmatrix}}_{\bm{\dot x}(t)} =\underbrace{\begin{bmatrix} 0 & 1 &0 \\ 0 & -\frac{K_f}{J}& \frac{K}{J} \\0 & -\frac{K}{L}& -\frac{R}{L} \end{bmatrix}}_{A} \underbrace{\begin{bmatrix} \theta(t) \\ \omega(t) \\ i(t)\end{bmatrix}}_{\bm{x}(t)} +\underbrace{\begin{bmatrix} 0 \\ 0 \\ \frac{1}{L}\end{bmatrix}}_{B}u(t)x˙(t)⎣⎡θ˙(t)ω˙(t)i˙(t)⎦⎤=A⎣⎡0001−JKf−LK0JK−LR⎦⎤x(t)⎣⎡θ(t)ω(t)i(t)⎦⎤+B⎣⎡00L1⎦⎤u(t)

输出方程为：

y(t)=[100]⎵C[θ(t)ω(t)i(t)]⎵x(t)y(t)=\underbrace{\begin{bmatrix} 1&0&0 \end{bmatrix}}_{C} \underbrace{\begin{bmatrix} \theta(t)\\ \omega(t)\\ i(t) \end{bmatrix}}_{\bm{x}(t)}y(t)=C[100]x(t)⎣⎡θ(t)ω(t)i(t)⎦⎤

在MATLAB里面构造状态方程也很简单，只要定义 A,B,C,DA , B , C , DA,B,C,D 四个矩阵就行了。

J = 3.2284E-6;%inetia
Kf = 3.5077E-6; %damping
K = 0.0274;  % torque constant
R = 4;   %resistance
L = 2.75E-6; %inductanceA = [0 1 00 -Kf/J K/J0 -K/L -R/L];
B = [0 ; 0 ; 1/L];
C = [1  0  0];
D = 0;
motor= ss(A,B,C,D);

我们也可以用simulink建模，具体如下：

当然，最简单的还是simulink的自带模块state space。

四、如何求解状态方程

方程我们有了，接下来就是如何求解了，我们先来审视一下状态方程：

x˙(t)=Ax(t)+Bu(t)\bm{\dot x}(t)=A\bm{x}(t)+B\bm{u}(t)x˙(t)=Ax(t)+Bu(t)

很遗憾，这个方程不是很好解，我们先来个简化版的，如果输入为零， u(t)=0\bm{u}(t) =0u(t)=0 ，会是怎么样呢？此时 x˙(t)=Ax(t)\bm{\dot x}(t)=A\bm{x}(t)x˙(t)=Ax(t) ，借鉴我们我们在第二节的思路，我们知道， x˙(t)=λx(t){\dot x}(t)=\lambda{x}(t)x˙(t)=λx(t) 的解为： x(t)=x(0)eλtx(t)=x(0)e^{\lambda t}x(t)=x(0)eλt ，类比一下， x˙(t)=Ax(t)\bm{\dot x}(t)=A\bm{x}(t)x˙(t)=Ax(t) 的解为：

x(t)=x(0)eAt\bm{x}(t)=\bm{x}(0)e^{At}x(t)=x(0)eAt

当然这可以通过严谨的数学求解得到。可见，在零输入条件下，系统的输出由初始条件 x(0)\bm{x}(0)x(0) 以及矩阵 AAA 决定。这个结论非常重要，是线性状态方程的核心概念，后续状态反馈以及观测器模型处理方式都是将有输入的方程转化成零输入的方程。

我们自来看一下 u(t)≠0\bm{u}(t) \ne 0u(t)̸=0 时状态方程的解，这个稍微复杂一点，解的过程不重要，我们直接给结论：

x(t)=eAtx(0)⎵Comp. func+eAtq(t)⎵Part. int\bm{x(t)}=\underbrace{e^{At}\bm{x}(0)}_{\text{Comp. func}}+\underbrace{ e^{At}q(t)}_{\text{Part. int}}x(t)=Comp. funceAtx(0)+Part. inteAtq(t)

第一项我们称之为通解，它是自由方程（无输入） x˙(t)=Ax(t)\bm{\dot x}(t)=A\bm{x}(t)x˙(t)=Ax(t) 的解，它取决于初始状态 x(0)\bm{x}(0)x(0) 以及系统的模态（ e−λ1t,e−λ2t,⋯e^{-\lambda_1 t},e^{-\lambda_2 t},\cdotse−λ1t,e−λ2t,⋯ 等，即A 对应的特征函数），这个好理解，我们前面说过了。第二项是特解，也就是系统初始条件为零时（x(0)=0\bm{x}(0)=0x(0)=0 ）剩余的那部分解，它取决于系统的输入 u(t)\bm{u}(t)u(t) 以及系统的模态 e−λ1t,e−λ2t,⋯e^{-\lambda_1 t},e^{-\lambda_2 t},\cdotse−λ1t,e−λ2t,⋯ 。

把 x(t)=eAtq(t)\bm{x}(t)=e^{At}\bm{q}(t)x(t)=eAtq(t) 代入状态方程 x˙(t)=Ax(t)+Bu(t)\bm{\dot x}(t)=A\bm{x}(t)+B\bm{u}(t)x˙(t)=Ax(t)+Bu(t) 可得：

AeAtq(t)+eAtq˙(t)=AeAtq(t)+Bu(t)Ae^{At}\bm{q}(t)+e^{At}\dot {\bm{q}}(t)=Ae^{At}\bm{q}(t)+B\bm{u}(t)AeAtq(t)+eAtq˙(t)=AeAtq(t)+Bu(t)

即 q˙(t)=e−AtBu(t)\dot {\bm{q}}(t)=e^{-At}B\bm{u}(t)q˙(t)=e−AtBu(t)

两边积分： q(t)=∫0te−AτBu(τ)dτ{\bm{q}}(t)=\int_0^te^{-A\tau}B\bm{u}(\tau)d\tauq(t)=∫0te−AτBu(τ)dτ

于是我们可以求得状态方程总的解为：

x(t)=eAtx(0)⎵Comp. func+∫0teA(t−τ)Bu(τ)dτ⎵Part. int\bm{x(t)}=\underbrace{e^{At}\bm{x}(0)}_{\text{Comp. func}}+\underbrace{\int_0^{t} e^{A(t-\tau)}B\bm{u}(\tau)d\tau}_{\text{Part. int}}x(t)=Comp. funceAtx(0)+Part. int∫0teA(t−τ)Bu(τ)dτ

同样的方法，不难求得输出方程的解为：

y(t)=CeAtx(0)⎵Comp. int+∫0tCeA(t−τ)Bu(τ)dτ⎵Part. int+Du(t)⎵Feedforward\bm{y(t)}=\underbrace{Ce^{At}\bm{x}(0)}_{\text{Comp. int}}+\underbrace{\int_0^{t} Ce^{A(t-\tau)}B\bm{u}(\tau)d\tau }_{\text{Part. int}}+\underbrace{D\bm{u}(t)}_{\text{Feedforward}}y(t)=Comp. intCeAtx(0)+Part. int∫0tCeA(t−τ)Bu(τ)dτ+FeedforwardDu(t)

可见，对于状态方程，无论是特解还是通解， eAte^{At}eAt 都是非常重要的组成部分，那这个函数又是什么样的呢？我们知道， ete^{t}et 的泰勒展开为：

et=1+t+t22!+t33!+t44!+⋯e^t=1+t+\frac{t^2}{2!}+\frac{t^3}{3!}+\frac{t^4}{4!}+\cdotset=1+t+2!t2+3!t3+4!t4+⋯

类比一下， eAte^{At}eAt 则可以展开为：

eAt=I+At+(At)22!+(At)33!+(At)44!+⋯e^{At}=I+At+\frac{{(At)}^2}{2!}+\frac{{(At)}^3}{3!}+\frac{{(At)}^4}{4!}+\cdotseAt=I+At+2!(At)2+3!(At)3+4!(At)4+⋯

又 v−1Av=λ，所以A=vλv−1\bm{v}^{-1}A\bm{v}=\bm{\lambda} ，所以 A=\bm{v}\bm{\lambda}\bm{v}^{-1}v−1Av=λ，所以A=vλv−1 ，代入上面泰勒展开：

eAt=v(I+λt+(λt)22!+(λt)33!+(λt)44!+⋯)v−1=veλtv−1\begin{aligned}e^{At}&=\bm{v}\Big(I+\bm{\lambda}t+\frac{{(\bm{\lambda}t)}^2}{2!}+\frac{{(\bm{\lambda}t)}^3}{3!}+\frac{{(\bm{\lambda}t)}^4}{4!}+\cdots\Big)\bm{v}^{-1}\\ &=\bm{v}e^{\bm{\lambda}t}\bm{v}^{-1} \end{aligned} eAt=v(I+λt+2!(λt)2+3!(λt)3+4!(λt)4+⋯)v−1=veλtv−1

可见，状态方程的解就是矩阵 AAA 特征值 λ\bm{\lambda}λ 对应的模态 eλte^{\bm{\lambda}t}eλt 的叠加，因此，想要系统的稳定，特征值 λ\bm{\lambda}λ 必须位于复平面的负半部分，否则，对应的模态将会发散。

举个简单的例子，假设 A=aA=aA=a ， B=1B=1B=1 ， C=1C=1C=1 ， D=0D=0D=0 ，这就是一个最简单一阶单输入单输出系统，假如这个系统的输如为阶跃信号，即

u(t)={0t<01t≥0u(t)=\left\{ \begin{aligned} &0 &t<0\\ &1 &t\geq0\\ \end{aligned} \right.u(t)={01t<0t≥0

系统的初始状态为零，即 x(0)=0x(0)=0x(0)=0 ，那么这个系统的响应是啥呢？直接代入上面的公式就行了：

y(t)=∫0tea(t−τ)dτ=1−eaty(t)=\int_0^t e^{a(t-\tau)}d\tau=1-e^{at}y(t)=∫0tea(t−τ)dτ=1−eat

非常简单，而且最重要的是，我们不需要进行拉普拉斯变换了，直接在时域就可以求解。因此，这种方法也称之为时域法。

很多时候，我们都是零初始条件，因此只剩特解为（单位阶跃响应）：

y(t)=∫0tCeA(t−τ)Bu(τ)dτ+Du(t)=∫0tCeA(t−τ)Bdτ+D\begin{aligned} \bm{y}(t)&=\int_0^tCe^{A(t-\tau)}B{u}(\tau)d\tau+D{u}(t)\\ &=\int_0^tCe^{A(t-\tau)}Bd\tau+D\\ \end{aligned}y(t)=∫0tCeA(t−τ)Bu(τ)dτ+Du(t)=∫0tCeA(t−τ)Bdτ+D

求解结果为：

y(t)=CA−1eAtB⎵Transient+D−CA−1B⎵Steady State\bm{y}(t)=\underbrace{CA^{-1}e^{At}B}_{\text{Transient}}+\underbrace{D-CA^{-1}B}_{\text{Steady State}}y(t)=TransientCA−1eAtB+Steady StateD−CA−1B

还是前面直流电机的例子，假设我们要观察速度输出，即将 C=[100]C= \begin{bmatrix} 1&0&0\end{bmatrix}C=[100] 改为 C=[010]C= \begin{bmatrix} 0&1&0\end{bmatrix}C=[010] ，我们直接在时域计算输出，就可以的到电机的在单位阶跃下的速度响应了。

由 y(t)\bm{y}(t)y(t) 解的方程可以看出，假定系统是稳定的，在零初始条件下，解分成两部分，一部分瞬态解，随时间变化，逐渐衰减到零，另外一部分为稳态解。显然，稳态解是我们想要的，而瞬态解，我们则希望它尽快衰减到零，这就要求系统的模态 e−λ1t,e−λ2t,⋯e^{-\lambda_1 t},e^{-\lambda_2 t},\cdotse−λ1t,e−λ2t,⋯ 要尽快衰减，这可以通过设定矩阵 AAA 的特征值来实现，那为什么很多书上都写要通过系统的极点配置来实现呢？——两者本质一回事，我们来推导一下。

我们先来看看状态方程和传递函数的关系。老办法，还是拉普拉斯变换。因为前馈比较容易分析，我们先忽略 Du(t)D\bm{u}(t)Du(t) 。对状态方程两边同时进行拉普拉斯变换：

sx(s)−x(0)=Ax(s)+Bu(s)s\bm{x}(s)-\bm{x}(0)=A\bm{x}(s)+B\bm{u}(s)sx(s)−x(0)=Ax(s)+Bu(s)

适当的变形可以得到：

x(s)=(sI−A)−1(x(0)+Bu(s))\bm{x}(s)=(sI-A)^{-1}\Big(\bm{x}(0)+B\bm{u}(s)\Big)x(s)=(sI−A)−1(x(0)+Bu(s))

由于我们计算传递函数时一般忽略初始条件的，所以：

x(s)=(sI−A)−1Bu(s)\bm{x}(s)=(sI-A)^{-1}B\bm{u}(s)x(s)=(sI−A)−1Bu(s)

对输出方程两边同时进行拉普拉斯变换：

y(s)=Cx(s)\bm{y}(s)=C\bm{x}(s)y(s)=Cx(s)

所以系统的传递函数为：

y(s)u(s)=C(sI−A)−1B=C(sI−A)∗B∣sI−A∣\frac{\bm{y}(s)}{\bm{u}(s)}=C(sI-A)^{-1}B=\frac{C(sI-A)^*B}{\left| sI-A \right|}u(s)y(s)=C(sI−A)−1B=∣sI−A∣C(sI−A)∗B

其中 (sI−A)∗为(sI−A)(sI-A)^* 为 (sI-A)(sI−A)∗为(sI−A) 伴随矩阵， ∣sI−A∣\left| sI-A \right|∣sI−A∣ 为 (sI−A)(sI-A)(sI−A) 行列式。而矩阵 A 的特征值怎么计算呢？

∣λI−A∣=0\left| \lambda I-A \right|=0∣λI−A∣=0

把特征值 λ\lambdaλ 换成 sss 两个样子就一样了，也就是说矩阵 AAA 的特征值是传递函数的极点。在小潘的另一篇文章J Pan：如何入门自动控制理论中，我们分析了，传递函数的极点，就代表了系统的一种模态（ eλte^{\lambda t}eλt ），而这与本文中所说的微分算子 ddt\frac{d}{dt}dtd 特征函数就是一回事，因此是统一的。因此，对于控制系统，合理的极点配置，对于系统的动态响应是非常重要的。

五、可控性和可观性

假设我们一个状态方程如下：

x˙(t)=[2−3−1]⎵Ax(t)+[001310]⎵B[u1(t)u2(t)]\bm{\dot x}(t)= \underbrace{\begin{bmatrix} 2 & &\\ &-3&\\ & &-1 \end{bmatrix}}_{\text{A}} \bm{x}(t) +\underbrace{\begin{bmatrix} 0&0\\ 1&3\\ 1&0 \end{bmatrix}}_{\text{B}}\begin{bmatrix} u_1(t)\\ u_2(t)\end{bmatrix}x˙(t)=A⎣⎡2−3−1⎦⎤x(t)+B⎣⎡011030⎦⎤[u1(t)u2(t)]

很显然，这个系统有3个模态将会在瞬态响应中出现，分别为 e2t,e−t和e−3te^{2t} ,e^{-t} 和e^{-3t}e2t,e−t和e−3t ，显然， e−te^{-t}e−t 和 e−3te^{-3t}e−3t 是衰减的模态，而 e2te^{2t}e2t则是发散的模态。由于矩阵 BBB 的第一行为 [00][0 \quad 0][00] ，因此这个模态是不受输入 u1(t)和u2(t)u_1(t) 和 u_2(t)u1(t)和u2(t) 控制的，可以想象的到，这个系统会很危险——发散的！

通过前面我们分析，直觉告诉我们，是不是有模态不可控，貌似和矩阵AAA （确定模态）即矩阵 BBB （确定输入）有关，那怎么才能通过这两个矩阵来判断系统是否可控呢？经过数学推导，我们发现可以这样进行简单的判断，构造矩阵：

C=[BABA2B⋯An−1B]C=[B\quad AB \quad A^2B \quad \cdots A^{n-1}B]C=[BABA2B⋯An−1B]

当这个矩阵的秩为 nnn 时，系统是可控的，否则不可控。当然，我们一般不用手算，MATLAB代码为：

Co = ctrb(A,B);
rank(Co)

我们再来看看可观性，这要从输出方程来看：

[y1(t)y2(t)]=[021010]⎵C[x1(t)x2(t)x3(t)]\begin{bmatrix} y_1(t) \\ y_2(t)\end{bmatrix} =\underbrace{\begin{bmatrix} 0&2&1\\ 0&1&0\\ \end{bmatrix}}_{\text{C}}\begin{bmatrix} x_1(t)\\ x_2(t)\\ x_3(t)\end{bmatrix}[y1(t)y2(t)]=C[002110]⎣⎡x1(t)x2(t)x3(t)⎦⎤

状态方程还和前面一致，系统有三个模态e−t,e−3te^{-t} , e^{-3t}e−t,e−3t 和 e2te^{2t}e2t ，因为矩阵 CCC 中第一列为零，也就是说 x1(t)x_1(t)x1(t) 对应的模态 e2te^{2t}e2t ,它既未在输出变量 y1(t)y_1(t)y1(t) 中出现，也未在 y2(t)y_2(t)y2(t) 中出现，因此状态变量x1(t)x_1(t)x1(t) 是不可观测的。

同样，系统的可观性也可通过矩阵的秩来判断：

O=[CACA2C⋮An−1C]O=\begin{bmatrix} C \\ AC \\ A^2C \\ \vdots \\ A^{n-1}C\end{bmatrix}O=⎣⎢⎢⎢⎢⎢⎡CACA2C⋮An−1C⎦⎥⎥⎥⎥⎥⎤

如果矩阵 OOO 的秩为 nnn ，则系统可观，否则不可观。MATLAB代码为：

O = obsv(A,C);
rank(O)

六、如何用状态反馈进行极点配置

前面我们在说到状态方程的时域解时说到了，无论是通解还是特解，都取决于 eAte^{At}eAt ，保证系统的稳定，就要求矩阵 AAA 的特征值也就是极点位于复平面的左半部分。当然，要想系统有更好的响应，需要对极点做更优化的配置，这就需要用到状态反馈了。

具体做法就是将系统的全状态引入到输入里面去，见下图。

其中，KcK_cKc 所在的支路就是提供状态反馈，由于该条之路的存在，状态方程发生了改变。此时

u(t)=−Kcx(t)+r(t)\bm{u}(t)=-\bm{K}_c\bm{x}(t)+ \bm{r}(t)u(t)=−Kcx(t)+r(t)

则状态方程变为：

x˙(t)=(A−BKc)⎵new Ax(t)+Br(t)\bm{\dot x}(t)=\underbrace{(A-BK_c)}_{\text{new A}} \bm{x}(t)+B\bm{r}(t)x˙(t)=new A(A−BKc)x(t)+Br(t)

因此，闭环控制系统的极点分布为：

∣λI−(A−BKc)∣=0\left| \lambda I-(A-BK_c) \right|=0∣λI−(A−BKc)∣=0

因此，我们想要把极点配置在任何位置，只需要改变矩阵 KcK_cKc 即可。

还是前面DC motor的例子，假设我们想把极点配置在 [−100±100i−200]\begin{bmatrix} -100\pm 100i &-200\end{bmatrix}[−100±100i−200] ，我们借助MATLAB函数place帮我们计算一下，具体代码如下：

p1 = -100+100i;
p2 = -100-100i;
p3 = -200;
K = place(A,B,[p1, p2, p3])
t = 0:0.001:0.05;
sys_cl = ss(A-B*K,B,C,D);
step(sys_cl,t)

此时，MATLAB给出的状态反馈增益为：Kc=[2.0000−0.0704−0.6020]K_c= \begin{bmatrix} 2.0000 &-0.0704 &-0.6020\end{bmatrix}Kc=[2.0000−0.0704−0.6020]

我们再搭个Simulink模型验证一下：

位置响应曲线如下。可见，瞬态部分最终完全衰减，最终只剩下稳态的部分。但注意：此时系统的输出和输入还有较大的差别，该如何解决？

可见，单位阶跃信号的响应样子和阶跃信号差不多了，但是幅值差了很多倍，我们需要通过一个增益调整过来，我们把框图改成如下形式，在指令通道上增加了 krk_rkr ，那该如何计算这个增益呢？

我们从稳态入手，当系统到达稳态的时候，会出现什么呢？——什么是稳态呢？就是状态量不再变化，换个说法，也就是状态量的导数为零。我们来看一下状态方程：

x˙(t)=(A−BKc)x(t)+Bkrr(t)\bm{\dot x}(t)=(A-BK_c) \bm{x}(t)+Bk_r\bm{r}(t)x˙(t)=(A−BKc)x(t)+Bkrr(t)

当 x˙e(t)=0\bm{\dot x_e}(t)=0x˙e(t)=0 时，我们把 x˙e(t)\bm{\dot x_e}(t)x˙e(t) 称之为平衡点，这时可以得到：

(A−BKc)xe(t)+Bkrr(t)=0(A-BK_c) \bm{x_e}(t)+Bk_r\bm{r}(t)=0(A−BKc)xe(t)+Bkrr(t)=0

所以

xe(t)=−(A−BKc)−1Bkrr(t)\bm{x_e}(t)=-(A-BK_c)^{-1}Bk_r\bm{r}(t)xe(t)=−(A−BKc)−1Bkrr(t)

如果我们想要输入输出一致，则有 r(t)=Cxe(t)\bm{r}(t)=C \bm{x_e}(t)r(t)=Cxe(t) ，代入上式就可以得到;

kr=−1/(C(A−BKc)−1B)k_r=-1/\big(C(A-BK_c)^{-1}B\big)kr=−1/(C(A−BKc)−1B)

我们搭个模型试一下：

输出曲线如下，可见，输出基本跟随输入，但是放大了看，还是有稳态误差。

七、如何进行跟踪控制

前面一节我们看到了，只进行状态反馈，只能调整系统模态的衰减快慢，虽然能保证输入输出一致，但仍存在在稳态误差，而且因为没有对输出进行反馈控制，系统的抗干扰能力很差，也就是说，还不能算是一个很好的控制系统。那怎么办呢？——在经典控制理论里面，我们是怎么消除静态误差的呢？——用积分，那有什么借鉴意义呢？——我们是不是可以增加一个新状态变量（state）来表征输出和输入之间误差？如果这个状态变量为零，那输入和输出就一致了。

状态变量定义如下：

xi(t)=∫y(t)−r(t)dt\bm{x_i}(t)=\int \bm{y}(t)-\bm{r}(t)dtxi(t)=∫y(t)−r(t)dt

两边都微分一下，就是：

x˙i(t)=y(t)−r(t)=Cx(t)−r(t)\begin{aligned}\bm{\dot x_i}(t)&= \bm{y}(t)-\bm{r}(t)\\&=C\bm{x}(t)-\bm{r}(t) \end{aligned}x˙i(t)=y(t)−r(t)=Cx(t)−r(t)

把这个新变量扩充到原来的状态方程

{x˙(t)=Ax(t)+Bu(t)y(t)=Cx(t)+Du(t)\left\{ \begin{array}{lcl} \bm{\dot x}(t)=A\bm{x}(t)+B\bm{u}(t) \\ \bm{y}(t)=C\bm{x}(t)+D\bm{u}(t) \end{array} \right.{x˙(t)=Ax(t)+Bu(t)y(t)=Cx(t)+Du(t)

可以得到：

[x˙(t)x˙i(t)]⎵x˙a(t)=[A0C0]⎵Aa[x(t)xi(t)]+[B0]⎵Bau(t)+[0−I]⎵Brr(t)\underbrace{ \begin{bmatrix} \bm{\dot x}(t) \\ \bm{\dot x_i}(t)\end{bmatrix} }_{\bm{\dot x}_a(t)} = \underbrace{ \begin{bmatrix} A &0\\C & 0\end{bmatrix} }_{A_a} \begin{bmatrix} \bm{ x}(t) \\ \bm{x_i}(t)\end{bmatrix} + \underbrace{\begin{bmatrix} B\\ 0\end{bmatrix}}_{{B_a}} \bm{u}(t) + \underbrace{\begin{bmatrix} 0\\ -I\end{bmatrix}}_{{B_r}} \bm{r}(t)x˙a(t)[x˙(t)x˙i(t)]=Aa[AC00][x(t)xi(t)]+Ba[B0]u(t)+Br[0−I]r(t)

我们把扩展的矩阵起个新的名字，分别为 Aa,Ba,BrA_a , B_a , B_rAa,Ba,Br 。还拿前面的DC motor的例子，我们把误差状态变量增加进去，就可以得到新的状态方程：

[θ˙(t)ω˙(t)i˙(t)ν˙(t)]=[01000−KfJKJ00−KL−RL01000][θ(t)ω(t)i(t)ν(t)]+[001L0]u(t)+[000−1]r(t)\begin{bmatrix} \dot\theta(t) \\ \dot\omega(t) \\ \dot i(t)\\\dot{\nu}(t)\end{bmatrix} =\begin{bmatrix} 0 & 1 &0 &0\\ 0 & -\frac{K_f}{J}& \frac{K}{J} &0\\0 & -\frac{K}{L}& -\frac{R}{L} &0\\1&0&0&0\end{bmatrix} \begin{bmatrix} \theta(t) \\ \omega(t) \\ i(t)\\\nu(t)\end{bmatrix} +\begin{bmatrix} 0 \\ 0 \\ \frac{1}{L}\\0\end{bmatrix}u(t) +\begin{bmatrix} 0 \\ 0 \\ 0\\-1\end{bmatrix}r(t)⎣⎢⎢⎡θ˙(t)ω˙(t)i˙(t)ν˙(t)⎦⎥⎥⎤=⎣⎢⎢⎡00011−JKf−LK00JK−LR00000⎦⎥⎥⎤⎣⎢⎢⎡θ(t)ω(t)i(t)ν(t)⎦⎥⎥⎤+⎣⎢⎢⎡00L10⎦⎥⎥⎤u(t)+⎣⎢⎢⎡000−1⎦⎥⎥⎤r(t)

以及新的输出方程：

y(t)=[1000]⎵Ca[θ(t)ω(t)i(t)ν(t)]y(t)=\underbrace{\begin{bmatrix} 1&0&0&0\end{bmatrix}}_{C_a} \begin{bmatrix} \theta(t) \\ \omega(t) \\ i(t)\\\nu(t)\end{bmatrix}y(t)=Ca[1000]⎣⎢⎢⎡θ(t)ω(t)i(t)ν(t)⎦⎥⎥⎤

式中 ν(t)\nu(t)ν(t) 就是前面的 xi(t)x_i(t)xi(t) 。

画成框图就容易多了：

接下来的任务就是该如何去确定反馈增益 KiK_iKi 。

根据上面的框图，可以得到输入 u(t)\bm{u}(t)u(t) 的表达式为：

u(t)=−Kcx(t)−Kiν(t)=−[KcKi][x(t)ν(t)]\bm{u}(t)=-K_c\bm{x}(t)-K_i\bm{\nu}(t) =-\begin{bmatrix} K_c &K_i\end{bmatrix} \begin{bmatrix} \bm{x}(t) \\\bm{\nu}(t)\end{bmatrix}u(t)=−Kcx(t)−Kiν(t)=−[KcKi][x(t)ν(t)]

令 xa(t)=[x(t)ν(t)]以及Ka=[KcKi]\bm{x_a}(t)= \begin{bmatrix} \bm{x}(t) \\\bm{\nu}(t)\end{bmatrix} 以及 K_a=\begin{bmatrix} K_c &K_i\end{bmatrix}xa(t)=[x(t)ν(t)]以及Ka=[KcKi]

则可以得到有输出反馈的系统状态方程：

x˙a(t)=(Aa−BaKa)⎵new Axa(t)+Brr(t)\bm{\dot x_a}(t)=\underbrace{(A_a-B_aK_a)}_{\text{new A}} \bm{x_a}(t)+B_r\bm{r}(t)x˙a(t)=new A(Aa−BaKa)xa(t)+Brr(t)

以及输出方程

y(t)=Caxa(t)\bm{y}(t)=C_a\bm{x_a}(t)y(t)=Caxa(t)

因此，闭环控制系统的极点分布为：

∣λI−(Aa−BaKa)∣=0\left| \lambda I-(A_a-B_aK_a) \right|=0∣λI−(Aa−BaKa)∣=0

也就是说，我们只要配置矩阵 Aa−BaKaA_a-B_aK_aAa−BaKa 的极点就可以了，具体代码如下：

Aa = [0 1 0 00 -Kf/J K/J 00 -K/L -R/L 01 0 0 0];
Ba = [0 ; 0 ; 1/L ; 0 ];
Br = [0 ; 0 ; 0; -1];
Ca = [1 0 0 0];
Da = [0];p4 = -300;
Ka = place(Aa,Ba,[p1,p2,p3,p4]);
Kc=Ka(1:3);
Ki=Ka(4);
t = 0:0.001:0.1;
sys_cl = ss(Aa-Ba*Ka,Br,Ca,Da);
step(sys_cl,t)

Simulink模型为：

系统的输出为：

可见，位置信号的动态响应能快速跟随单位阶跃信号，实现了很好的闭环控制。

注意，可能细心的童鞋可能注意到了，用积分的方式增加一个状态变量，此时系统的等效状态方程为：

x˙a(t)=(Aa−BaKa)⎵new Axa(t)+Br⎵new Br(t)\bm{\dot x_a}(t)=\underbrace{(A_a-B_aK_a)}_{\text{new A}} \bm{x_a}(t)+\underbrace{B_r}_{\text{new B}} \bm{r}(t)x˙a(t)=new A(Aa−BaKa)xa(t)+new BBrr(t)

此时等效系统并不是我们熟悉的零输入系统，也就是说，系统稳定时，状态变量 xa(t)\bm{x_a}(t)xa(t) 并不会衰减到零——这个没有关系，因为当输入 r(t)\bm{r}(t)r(t) 是阶跃信号时，如果 (Aa−BaKaA_a-B_aK_aAa−BaKa) 的极点都在复平面左半边，那 xa(t)\bm{x_a}(t)xa(t) 最终会稳定到一个恒定值，也就是 xi(t)=constant\bm{x_i}(t)=constantxi(t)=constant ，所以 x˙i(t)=y(t)−r(t)=0\bm{\dot x_i}(t)= \bm{y}(t)-\bm{r}(t)=0x˙i(t)=y(t)−r(t)=0 ，也就是 y(t)=r(t)\bm{y}(t)=\bm{r}(t)y(t)=r(t) ，因此，能够实现伺服跟踪控制。

注：这种处理方式只能保证被控对象Type 0型系统时阶跃激励（ 1/s1/s1/s ）稳态误差为零，至于Type 1及以上系统或者激励为斜坡（ 1/s21/s^21/s2 ）等系统，反馈部分的结构设计略有不同，感兴趣的请阅读相关专著。

八、全状态观测器

通过前面的介绍我们已经知道，状态反馈是一个强大的工具，能够实现任何预定的极点配置。但是如果仔细观察反馈状态变量

u(t)=−Kcx(t)\bm{u}(t)=-\bm{K}_c\bm{x}(t)u(t)=−Kcx(t)

就会发现仍有问题，那就是不是所有的变量都是可测量的，比如：

如果状态变量过多，比如有 x,x˙,x¨,θ,θ˙,θ¨,etcx,\dot x,\ddot x,\theta,\dot \theta,\ddot \theta,etcx,x˙,x¨,θ,θ˙,θ¨,etc ，这会需要大量的传感器；
很多时候 x˙,x¨\dot x,\ddot xx˙,x¨ 是从 xxx 微分得到的，会有很大的噪音；
有些工况无法使用传感器，比如说核反应堆。
那怎么办呢？——用观测器，那怎么样观测呢？一个很自然的想法就是我把原模型copy一份啊，来个一模一样的计算模型，这样各个变量不都可以算出来了吗？

比如我们可以把状态方程和输出方程复制一份

{x˙(t)=Ax(t)+Bu(t)y(t)=Cx(t)+Du(t)\left\{ \begin{array}{lcl} \bm{\dot x}(t)=A\bm{x}(t)+B\bm{u}(t) \\ \bm{y}(t)=C\bm{x}(t)+D\bm{u}(t) \end{array} \right.{x˙(t)=Ax(t)+Bu(t)y(t)=Cx(t)+Du(t)

然后就可以得到如下的示意图：

从框图可以看出，假如我们给两个模型（物理模型和数学模型）同样的输入，如果初始状态也相同，理论上两个模型的状态变量和输出应该是完全一样的。

实际上不可能这么完美，因为我们不可能知道完整的初始条件是什么样，因此不可能完全一样，而且两个模型也不可能完全一样，那输出有可能会差别很大。

控制的核心思想就是反馈，因此，当开环的模型存在较大误差，我们可以尝试一下是不是可以闭环控制，比如采用如下的形式：

这个框图啥意思呢？——就是将两个模型的输出量的误差反馈到观测模型的输入，这样只要输出有差别，观测器的状态就会进行跟随调整，保证两个模型的状态一致，进而保证输出也一致。这样观测器的方程就变为：

{x˙^(t)=Ax^(t)+Bu(t)+L(y(t)−y^(t))y^(t)=Cx^(t)+Du(t)\left\{ \begin{array}{lcl} \hat{\bm{\dot x}}(t)=A\hat{\bm{x}}(t)+B\bm{u}(t) +L\Big(\bm{y}(t)-\hat{\bm{y}}(t)\Big)\\ \hat{\bm{y}}(t)=C\hat{\bm{x}}(t)+D\bm{u}(t) \end{array} \right.{x˙^(t)=Ax^(t)+Bu(t)+L(y(t)−y^(t))y^(t)=Cx^(t)+Du(t)

将两个模型的输出方程带入到状态方程，就可以得到：

x˙^(t)=Ax^(t)+Bu(t)+LC(x(t)−x^(t))\hat{\bm{\dot x}}(t)=A\hat{\bm{x}}(t)+B\bm{u}(t) +LC\Big(\bm{x}(t)-\hat{\bm{x}}(t)\Big)x˙^(t)=Ax^(t)+Bu(t)+LC(x(t)−x^(t))

将观测器的状态方程与物理模型的状态方程相减：

x˙(t)−x˙^(t)=(A−LC)(x(t)−x^(t))\bm{\dot x}(t)-\hat{\bm{\dot x}}(t)=(A-LC)\Big({\bm{x}}(t)-\hat{\bm{x}}(t)\Big)x˙(t)−x˙^(t)=(A−LC)(x(t)−x^(t))

即

e˙(t)=(A−LC)e(t)\bm{\dot e}(t)=(A-LC)\bm{e}(t)e˙(t)=(A−LC)e(t)

所以，只要合理配置 (A−LC)(A-LC)(A−LC) 的极点即可。

还是之前的DC motor模型，因为观测器的数值作为反馈信号，这就要求观测的响应更快，为此，我们一般将 (A−LC)(A-LC)(A−LC) 的极点取到 (A−BKc)(A-BK_c)(A−BKc) 极点的10倍左右，具体MATLAB代码如下：

p1 = -100+100i;
p2 = -100-100i;
p3 = -200;
L = place(A',C',[p1*10, p2*10, p3*10])

完整的simulink模型为（考虑到实际情况，我们取 A^=0.9A\hat A=0.9AA^=0.9A）：

各状态变量及其观测值的响应为：

可见，观测值在很快的时间内收敛到真实值。整个系统的阶跃响应如下：

九、缩减状态观测器

无论需不需要，全状态观测器是对所有的状态进行观测。但是在有些系统中，很多状态是可以测量的，这是时候就有两个选择，要么用全状态观测器，然后选择用测量值或者观测值或者两者都用。另一个选择。使用缩减观测器，只对不能进行测量的状态进行观测。这部分公式比较多，但是推导过程展现的处理方式有很大的借鉴意义，有兴趣的建议按照思路自己推演一遍。

我们先要将状态分为可测量 xm(t)\bm{x_m}(t)xm(t) 和不可测量 xr(t)\bm{x_r}(t)xr(t) 两部分，然后用分块矩阵将状态方程改为：

[x˙m(t)x˙r(t)]=[A11A12A21A22][xm(t)xi(t)]+[B1B2]u(t)(1)\begin{bmatrix} \bm{\dot x_m}(t) \\ \bm{\dot x_r}(t)\end{bmatrix} =\begin{bmatrix} A_{11} &A_{12}\\A_{21} & A_{22}\end{bmatrix} \begin{bmatrix} \bm{ x_m}(t) \\ \bm{x_i}(t)\end{bmatrix} + \begin{bmatrix} B_1\\ B_2\end{bmatrix}\bm{u}(t) (1)[x˙m(t)x˙r(t)]=[A11A21A12A22][xm(t)xi(t)]+[B1B2]u(t)(1)

输出方程为：

y(t)=[C10][xm(t)xr(t)]\bm{y}(t)=\begin{bmatrix} C_1 & 0\end{bmatrix} \begin{bmatrix} \bm{x_m}(t) \\ \bm{x_r}(t)\end{bmatrix}y(t)=[C10][xm(t)xr(t)]

上一节我们说了，全状态观测器的方程为：

{x˙^(t)=Ax^(t)+Bu(t)+L(y(t)−y^(t))y^(t)=Cx^(t)+Du(t)(2)\left\{ \begin{array}{lcl} \hat{\bm{\dot x}}(t)=A\hat{\bm{x}}(t)+B\bm{u}(t) +L\Big(\bm{y}(t)-\hat{\bm{y}}(t)\Big)\\ \hat{\bm{y}}(t)=C\hat{\bm{x}}(t)+D\bm{u}(t) \end{array} \right. (2){x˙^(t)=Ax^(t)+Bu(t)+L(y(t)−y^(t))y^(t)=Cx^(t)+Du(t)(2)

缩减状态观测器处理方式基本一致，就是要找到自己的等效表达式，假设 D=0D=0D=0 ，缩减观测器的等效形式为：

{x˙r^(t)=Arxr^(t)+ur(t)+L(yr(t)−yr^(t))yr^(t)=Crxr^(t)\left\{ \begin{array}{lcl} \hat{\bm{\dot x_r}}(t)=A_r\hat{\bm{x_r}}(t)+\bm{u_r}(t) +L\Big(\bm{y_r}(t)-\hat{\bm{y_r}}(t)\Big)\\ \hat{\bm{y_r}}(t)=C_r\hat{\bm{x_r}}(t) \end{array} \right.{x˙r^(t)=Arxr^(t)+ur(t)+L(yr(t)−yr^(t))yr^(t)=Crxr^(t)
将状态方程（1）下半部分展开，得到：

x˙r(t)=A22⎵Arxr(t)+(A21xm(t)+B2u(t))⎵ur(t)(3)\begin{aligned} \bm{\dot x_r}(t)&=\underbrace{A_{22}}_{A_r}\bm{x_r}(t)+\underbrace{\Big(A_{21}\bm{x_m}(t)+B_2\bm{u}(t)\Big)}_{\bm{u_r}(t)}\\ \end{aligned} (3)x˙r(t)=ArA22xr(t)+ur(t)(A21xm(t)+B2u(t))(3)

可见，这和我们想要的缩减状态观测器基本一致，于是我们得到缩减状态观测器的 ArA_rAr 和等效输出ur(t)\bm{u_r}(t)ur(t) 。再将状态方程（1）的上半部分展开，得到：

x˙m(t)=A11xm(t)+A12xr(t)+B1u(t)\begin{aligned} \bm{\dot x_m}(t)&=A_{11}\bm{x_m}(t)+A_{12}\bm{x_r}(t)+B_1\bm{u}(t)\\ \end{aligned}x˙m(t)=A11xm(t)+A12xr(t)+B1u(t)

我们稍微改写一下形式：

A12xr(t)=x˙m(t)−A11xm(t)−B1u(t)⎵yr(t)A_{12}\bm{x_r}(t)=\underbrace{\bm{\dot x_m}(t)-A_{11}\bm{x_m}(t)-B_1\bm{u}(t)}_{\bm{y_r}(t)}A12xr(t)=yr(t)x˙m(t)−A11xm(t)−B1u(t)

上式右半部分取决于状态变量 xr(t)\bm{x_r}(t)xr(t) ，因此可以看成是输出 yr(t){\bm{y_r}(t)}yr(t) ，这样我们就有了输出方程。因此，我们就可以得到缩减状态观测器：

x˙^r(t)=A22x^r(t)+ur(t)+L(yr(t)−y^r(t))（4）{\bm{\hat {\dot x}_r}}(t)=A_{22}{\bm{\hat {x}_r}}(t)+\bm{u_r}(t)+L\Big(\bm{y_r}(t)-\bm{\hat {y}_r}(t)\Big) （4）x˙^r(t)=A22x^r(t)+ur(t)+L(yr(t)−y^r(t))（4）

将 yr(t)=x˙m(t)−A11xm(t)−B1u(t)\bm{y_r}(t)=\bm{\dot x_m}(t)-A_{11}\bm{x_m}(t)-B_1\bm{u}(t)yr(t)=x˙m(t)−A11xm(t)−B1u(t)

以及 y^r(t)=A12x^r(t)\bm{\hat {y}_r}(t)=A_{12} {\bm{\hat x}_r}(t)y^r(t)=A12x^r(t)

代入缩减状态观测器可以得到：

x˙^r(t)=A22x^r(t)+(A21xm(t)+B2u(t))⎵ur(t)+L((x˙m(t)−A11xm(t)−B1u(t)⎵yr(t)−A12x^r(t)⎵y^r(t))\begin{aligned} {\bm{\hat {\dot x}_r}}(t)&=A_{22}{\bm{\hat {x}_r}}(t)+\underbrace{\Big(A_{21}\bm{x_m}(t)+B_2\bm{u}(t)\Big)}_{{\bm{u_r}(t)}}\\&+L\underbrace{\Big((\bm{\dot x_m}(t)-A_{11}\bm{x_m}(t)-B_1\bm{u}(t)}_{\bm{y_r}(t)}-\underbrace{A_{12}\bm{\hat {x}_r}(t)}_{\bm{\hat {y}_r}(t)}\Big)\end{aligned}x˙^r(t)=A22x^r(t)+ur(t)(A21xm(t)+B2u(t))+Lyr(t)((x˙m(t)−A11xm(t)−B1u(t)−y^r(t)A12x^r(t))

重新整理一下，合并同类项，得到：

x˙^r(t)=(A22−LA12)x^r(t)+(A21−LA11)xm(t)⋯+(B2−LB1)u(t)+Lx˙m(t)\begin{aligned}\bm{\hat {\dot x}_r}(t) &=(A_{22}-LA_{12})\bm {\hat{x}_r}(t)+(A_{21}-LA_{11})\bm {{x}}_m(t) \cdots\\ \\&+(B_2-LB_1)\bm{u}(t)+L\bm{{\dot x}_m}(t) \end{aligned}x˙^r(t)=(A22−LA12)x^r(t)+(A21−LA11)xm(t)⋯+(B2−LB1)u(t)+Lx˙m(t)

仔细观察上式，我们发现 xm(t)\bm {{x}}_m(t)xm(t) 和 u(t)\bm{u}(t)u(t) 都是可测量或已知的，都没问题，但是 x˙m(t)\bm{{\dot x}_m}(t)x˙m(t) 未测量啊，怎么办？我们不妨把它归到xr(t)\bm {{x}}_r(t)xr(t) 里面，也就是说可以把这两个变量合成一个新变量：

x˙c(t)=x˙^r(t)−Lx˙m(t)\bm{\dot x_c}(t)=\bm{\hat {\dot x}_r}(t)-L\bm{{\dot x}_m}(t)x˙c(t)=x˙^r(t)−Lx˙m(t)

xc(t)=x^r(t)−Lxm(t)\bm{ x_c}(t)=\bm{\hat {x}_r}(t)-L\bm{{x}_m}(t)xc(t)=x^r(t)−Lxm(t)

这样，缩减状态观测器就变为：

x˙c(t)=(A22−LA12)x^r(t)+(A21−LA11)x^m(t)⋯+(B2−LB1)u(t)\begin{aligned} \bm{\dot x_c}(t)&=(A_{22}-LA_{12})\bm {\hat{x}_r}(t)+(A_{21}-LA_{11})\bm {\hat{x}}_m(t)\cdots\\ \\&+(B_2-LB_1)\bm{u}(t) \end{aligned}x˙c(t)=(A22−LA12)x^r(t)+(A21−LA11)x^m(t)⋯+(B2−LB1)u(t)

控制框图如下：

下一个问题是，我们该如何设计状态反馈矩阵 LLL ？

yr(t)=A12xr(t)\bm{y_r}(t)=A_{12}\bm{x_r}(t)yr(t)=A12xr(t)

代入缩减转台估计方程（4）可以得到：

x˙^r(t)=A22x^r(t)+(A21xm(t)+B2u(t))⎵ur(t)+LA12((xr(t)−x^r(t))⎵Error{\bm{\hat {\dot x}_r}}(t)=A_{22}{\bm{\hat {x}_r}}(t)+\underbrace{\Big(A_{21}\bm{x_m}(t)+B_2\bm{u}(t)\Big)}_{\bm{u_r}(t)}+LA_{12}\underbrace{\Big((\bm{x_r}(t)-\bm{\hat {x}_r}(t)\Big)}_{\text{Error}}x˙^r(t)=A22x^r(t)+ur(t)(A21xm(t)+B2u(t))+LA12Error((xr(t)−x^r(t))

与状态方程（3）相减，可以得到;
x˙r(t)−x˙^r(t)=(A22−LA12)(xr(t)−x^r(t)){\bm{{\dot x}_r}}(t)-{\bm{\hat {\dot x}_r}}(t)=(A_{22}-LA_{12})\Big ({\bm{{x}_r}}(t)-{\bm{\hat {x}_r}}(t)\Big)x˙r(t)−x˙^r(t)=(A22−LA12)(xr(t)−x^r(t))

即 e˙r(t)=(A22−LA12)er(t){\bm{{\dot e}_r}}(t)=(A_{22}-LA_{12}){\bm{{e}_r}}(t)e˙r(t)=(A22−LA12)er(t)

因此，只需要按照要求配置矩阵 (A22−LA12)(A_{22}-LA_{12})(A22−LA12) 的极点即可。

还是以之前说的DC motor为例，假设 xm(t)=[θ(t)ω(t)]\bm{x_m}(t)=\begin{bmatrix} \theta(t) \\ \omega(t)\end{bmatrix}xm(t)=[θ(t)ω(t)] ，即位置和速度是可观测的， xr(t)=i(t)\bm{x_r}(t)=i(t)xr(t)=i(t) ，电流是需要观测的。同样假设观测器模型是有误差的， A11,A12,A21,A22A_{11},A_{12},A_{21},A_{22}A11,A12,A21,A22 的具体数值如如下代码：

A11=A(1:2,1:2)*0.9;
A12=A(1:2,3)*0.9;
A21=A(3,1:2)*0.95;
A22=A(3,3)*0.9;
p5 = -200;
L = place(A22',A12',[p5*10])'

搭建完整的simulink模型为：

电流的观测值和实际值曲线见下图，可见即使观测器模型我们设置了较大的误差，电流仍得到了较好的估计。

整个系统的阶跃响应为：

十、小结
正所谓横看成岭侧成峰，无论是经典控制理论，还是现代控制理论，都是用来求解微分方程，因此在本质上应该是相通的。现代控制理论提供了一个非常强大的工具去解决控制问题，已经成为现代控制领域的理论基石，但是大多数教材过分强调其计算过程，现代控制理论背后的思想反而湮没在矩阵来回的变换中，最终导致大家始终不得要领。写这个文章的目的主要是梳理一下自己学习这们课程的思路，如果能对大家有一点帮助，那就善莫大焉了。当然，这一篇文章是不可能说明白所有的事，还有很多非常重要的内容，可能会需要一些的数学基础，比如稳定性问题，最优控制，kalman滤波，这些都需要在熟悉这篇文章的基础上才能有较好的体会，如果大家感兴趣，后续还可以再整理一些。

如何入门现代控制理论相关推荐

【控制】如何入门现代控制理论
在知乎看到的神文章,将自控和现控完全的串了一篇,太好的文章了,必须搬运分享过来,也感谢作者如何入门现代控制理论如何入门现代控制理论一.什么是特征值和特征向量二.什么是特征值和特征函数三.什么 ...
【控制】如何入门自动控制理论
在知乎看到的神文章,将自控和现控完全的串了一篇,太好的文章了,必须搬运分享过来,也感谢作者如何入门自动控制理论如何入门自动控制理论一.一个典型的系统响应有什么特点二.幅频响应和相频响应三.传 ...
如何入门自动控制理论
一.一个典型的系统响应有什么特点二.幅频响应和相频响应三.传递函数的零.极点是怎么个回事四.什么是闭环控制系统 https://zhuanlan.zhihu.com/p/42615269
c++矩阵作为函数输入变量_现代控制理论线性系统入门(七)输入输出解耦的控制器设计...
上一章传送门: 善道:线性系统控制入门(六)用能控标准型设计控制器zhuanlan.zhihu.com 在利用状态方程设计MIMO的能控标准型时,闭环系统的动态是完全不考虑输出变量而直接预给的,M ...
现代控制理论输出y_现代控制理论线性系统入门(三)输入输出变量的稳定性
上一篇传送门: 善道:线性系统控制入门(二)状态变量的解和渐近稳定性zhuanlan.zhihu.com 渐近稳定是状态变量的特性,如果再度考虑系统的输入输出行为,这个输入输出变量的稳定性可以由单变 ...
信号与线性系统管致中第六版pdf_现代控制理论线性系统入门(十)轨迹追踪反馈控制...
上一章传送门: 善道:线性系统控制入门(九)设计状态观测器zhuanlan.zhihu.com 上一章我们讨论了状态观测器的设计,实际上这是在实际操作中很必要的一个设计.因为实际上如果各个状态变量未 ...
simulink降维观测器设计_现代控制理论线性系统入门(九)设计状态观测器
上一章传送门: 善道:线性系统控制入门(八)最优控制--LQR控制器zhuanlan.zhihu.com 到目前为止,对于状态反馈控制器设计方法都是基于状态变量全是已知的.然而在很多情况下,这是不可 ...
职业：图像处理入门教程
公司里有哥们想学图像处理,认为这东西看起来挺简单,应该可以很快入门进阶.他人之心,吾忖度之,不能再忖度了,否则会变成以小人之心度君子之腹.所以,应哥们的要求,我给他介绍一个图像处理入门教程.看起来入门 ...
如何零基础入门FPGA？这篇文章让你吃透！
01 FPGA学习重点 1. 看代码,建模型只有在脑海中建立了一个个逻辑模型,理解FPGA内部逻辑结构实现的基础,才能明白为什么写Verilog和写C整体思路是不一样的,才能理解顺序执行语言和并行 ...

如何入门现代控制理论