C4 基于观测器的状态的状态反馈控制 Commande par retour d’état avec observateur

PLAN


  • 状态重构 Reconstruction d’état
    • 确定性情形:通过极点配置的观测器 Casdéterministe : Observateur par placement de pôles
    • 状态反馈控制和通过极点配置的观测器 Commandeparretourd’étatet observateur par placement de pôles
  • 优化方法 Approche optimale
    1. 卡尔曼滤波器(随机情形)Filtre de Kalman (cas stochastique)
    2. 线性二次型高斯(LQG)控制 CommandeLinéaire Quadratique Gaussienne (LQG)
  • 考虑未测量的恒定扰动 Prise encompted’une perturbation constante non mesurée

C4.1 状态重构 Reconstruction de l’état


状态重构指给定一个完全可观测的动态系统,使用观测器observateur状态重构器reconstructeur,能够在每一时刻根据系统的测量量(通常是输入和输出)来估计状态向量。

C4.1.1 观察器和观察系统


在很多情况下,\(x(t)\)无法测量,但\(y(t)\)可以。将输出\(y(t)\)视为传感器可测量的状态变量的集合,提出观测系统的动态结构如下,称\(\hat x\)为估计状态état estimé

\[ \begin{aligned}&\dot{\hat{x}}(t)=A \hat{x}(t)+B u(t)+L(y(t)-\hat{y}(t)) = (A-L C) \hat{x}(t)+B u(t)+L y(t)\\&\hat{y}(t)=C \hat{x}(t)\end{aligned} \]

该系统的输入为原系统输出\(y(t)\)和原系统输入\(u(t)\)。理想的观测系统应该具有合适的增益矩阵matrice de gains \(L \in \mathbb{R}^{n \times p}\),使得估计误差erreur d’estimation \(\tilde{x}(t)=x(t)-\hat{x}(t)\)趋近于0。


考虑误差\(\tilde x\),其演化方程为:

\[ \dot{\tilde{x}}(t)=(A-L C) \tilde{x}(t) \]

遵循自身的动态行为,而与状态的演化无关。因此,其稳定性取决于自身,等同于观测器的稳定性。回忆第二行关于稳定性的内容:

🪺 定理 1

线性时不变系统是稳定的当且仅当:

  • 矩阵A的所有特征值实部都非正:\(\forall \lambda \in \operatorname{sp}(A), \quad \operatorname{Re}(\lambda) \leq 0\)
  • 对于所有实部为零的特征值,其代数重数几何重数必须相等

线性时不变系统是渐进/指数稳定的,当且仅当:

  • A是Hurwitz 矩阵。

因此,需将系统的观测器通过 \((A, B, C)\) 的状态表示定义,并通过选择合适的 \(L\),使得矩阵 \(A-L C\) 的特征值为负数。此外,观测器的速度也十分重要。通常情况下,可以选择观测器的动态速度为被观测闭环系统速度的三倍。

C4.1.2 Wonham 观测定理


\(A \in \mathbb{R}^{n \times n} , C \in \mathbb{R}^{p \times n}\) 。对于任意期望的特征值集合 \(\left\{\lambda_1, \ldots, \lambda_n\right\}\) ,其中每个 \(\lambda_i\) 可以是真实数或复数(复数与其共轭成对出现),存在一个矩阵 \(L \in \mathbb{R}^{n \times p}\) ,使得:

$$

(A-L C)={_1, , _n}

$$

关于这一点,实际上可以基于第三章的Wonham定理结合可控和可测的对偶性证明。

CE4.1 观测伴随形式的例子


🪺 任何单变量可观测系统都可以用伴随形式表示,设演化矩阵和观测矩阵为:

\[ A=\left(\begin{array}{ccccc} 0 & 0 & \cdots & 0 & -a_0 \\ 1 & 0 & \cdots & 0 & -a_1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 0 & -a_{n-2} \\ 0 & 0 & \cdots & 1 & -a_{n-1} \end{array}\right), \quad C=\left(\begin{array}{lllll} 0 & 0 & \cdots & 0 & 1 \end{array}\right) \]

由于是单变量系统,增益矩阵\(L=\left(\ell_1, \ell_2, \ldots, \ell_n\right)^{\top} \in \mathbb{R}^n\),有:

\[ A-L C=\left(\begin{array}{ccccc} 0 & 0 & \cdots & 0 & -\left(a_0+\ell_1\right) \\ 1 & 0 & \cdots & 0 & -\left(a_1+\ell_2\right) \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 0 & -\left(a_{n-2}+\ell_{n-1}\right) \\ 0 & 0 & \cdots & 1 & -\left(a_{n-1}+\ell_n\right), \end{array}\right) \]

其特征多项式为:

\[ \operatorname{det}(\lambda I-(A-L C))=\lambda^n+\left(a_{n-1}+\ell_n\right) \lambda^{n-1}+\ldots+\left(a_1+\ell_2\right) \lambda+a_0+\ell_1 \]

对于任意一组特征值\(\left\{\lambda_i\right\}_{i \in\{1, \ldots, n\}}\)(实数或共轭复数对),存在唯一一组实系数\(\left\{p_i\right\}_{i \in\{1, \ldots, n\}}\),使得:

\[ \prod_{i=1}^n\left(\lambda-\lambda_i\right)=\lambda^n+\sum_{i=1}^n p_i \lambda^{\lambda^{i-1}} \]

由此可选择\(L\)中的元素。


C4.1.3 估计状态反馈控制 Commande par retour d'état estimé


对于状态:

\[ \begin{aligned} & \dot{x}(t)=A x(t)+B u(t) \\ & y(t)=C x(t) \end{aligned} \]

以及观测器:

\[ \dot{\hat{x}}(t)=A \hat{x}(t)+B u(t)+L(y(t)-C \hat{x}(t)) \]

在第三章中,我们设计闭环控制,使得\(u(t)=-K x(t)+My_c(t)\)。这建立在\(x(t)\)可以被获知的前提下。在本章中,使用估计状态替代状态向量,假设:

\[ u(t)=-K \hat{x}(t)+M y_c(t) \]

整个闭环系统,包括观测器,的状态方程可以被写作:

\[ \begin{aligned} {\left[\begin{array}{c} \dot{x}(t) \\ \dot{\hat{x}}(t) \end{array}\right] } & =\left[\begin{array}{c|c} A & -B K \\ \hline L C & A-B K-L C \end{array}\right]\left[\frac{x(t)}{\hat{x}(t)}\right]+\left[\frac{B M}{B M}\right] y_c(t) \\ y(t) & =\left[\begin{array}{l|l} C & 0 \end{array}\right]\left[\frac{x(t)}{\hat{x}(t)}\right] \end{aligned} \]


为了展示估计误差\(\tilde{x}(t)=x(t)-\hat{x}(t)\)的行为,使用一个新基重新表示系统:

\[ \begin{aligned} & {\left[\begin{array}{c} \dot{x}(t) \\ \hline \dot{\tilde{x}}(t) \end{array}\right]=\left[\begin{array}{c|c} A-B K & B K \\ \hline 0 & A-L C \end{array}\right]\left[\begin{array}{c} x(t) \\ \hline \tilde{x}(t) \end{array}\right]+\left[\begin{array}{c} B M \\ \hline 0 \end{array}\right] y_c(t)} \\ & y(t)=\left[\begin{array}{l|l} C & 0 \end{array}\right]\left[\begin{array}{l} x(t) \\ \tilde{x}(t) \end{array}\right] . \end{aligned} \]

这个公式表明,估计误差的动态行为与系统状态的演变无关,只取决于观测器的状态。

进一步观察,可以得到以下几点结论:

  • 系统的稳定性取决于全局演化矩阵\(\left[\begin{array}{c|c}A-B K & B K \\\hline 0 & A-L C\end{array}\right]\),而这是一个上三角分块矩阵,因此系统稳定性取决于标准闭环控制的演化矩阵\(A-BK\)的特征值和观察器演化矩阵\(A-LC\)的特征值。
  • 误差不可通过\(y_c(t)\)控制,自行演化。如果误差能初始化为\(0\),则始终为\(0\).
  • 基于估计状态的闭环控制系统的输入输出关系与标准的控制系统的输入输出关系相同。

分离原理

基于估计状态反馈的控制生成一个闭环系统,其特征值为未重构闭环系统的特征值,补充观察器的特征值。估计误差 \(\tilde{x}=x-\hat{x}\) 是自主的,且特别是不可通过输入 \(y_c\) 控制。闭环传递函数的输入 \(y_c\) 和输出之间的关系在有无估计的情况下完全相同。

值得注意的是,控制和观测是相互独立调整的,并且在施加基于观察器状态的控制时,它们不会改变。因此,基于估计状态反馈的控制设计可以分两步完成:首先调整基于标准闭环控制的控制,然后设计一个观察器,以确保状态估计足够快。

CE4.2 双积分器的例子


🪺 考虑双积分器的例子

其状态为:

\[ \binom{\dot{x}_1(t)}{\dot{x}_2(t)}=\left(\begin{array}{ll}0 & 1 \\0 & 0\end{array}\right)\binom{x_1(t)}{x_2(t)}+\binom{0}{1} u(t) \quad \text { et } \quad y(t)=x_1(t) \]

对于观察器,有:

\[ \begin{aligned}\dot{\hat{x}}(t)&=A \hat{x}(t)+B u(t)+L(y(t)-\hat{y}(t)) (t) \\\Rightarrow \binom{\dot{\hat{x}}_1(t)}{\hat{\hat{x}}_2(t)}&=\left(\begin{array}{ll} 0 & 1 \\ 0 & 0 \end{array}\right)\binom{\hat{x}_1(t)}{\hat{x}_2(t)}+\binom{0}{1} u(t)+\binom{l_1}{l_2}\left(y(t)-\hat{x}_1(t)\right)\\&=\left(\begin{array}{ll} -l_1 & 1 \\ -l_2 & 0 \end{array}\right)\binom{\hat{x}_1(t)}{\hat{x}_2(t)}+\binom{0}{1} u(t)+\binom{l_1}{l_2} y(t)\end{aligned} \]

对于二阶控制系统,有:

\[ \left|\begin{array}{cc} \lambda+l_1 & -1 \\ I_2 & \lambda \end{array}\right|=\lambda^2+l_1 \lambda+I_2 \equiv \lambda^2+2 \xi \omega_0 \lambda+\omega_0{ }^2 \rightarrow\left\{\begin{array}{c} l_1=2 \xi \omega_0 \\ l_2=\omega_0{ }^2 \end{array}\right. \]

C4.1.4 基于观察器的闭环系统的扰动抑制


可测量的情况

假设系统收到了扰动:

\[ \begin{gathered} \dot{x}(t)=A x(t)+B u(t)+B^{\prime} d(t) \\ y(t)=C x(t)+D^{\prime} d(t) \end{gathered} \]

此时,控制律:

\[ u(t)=-K x(t)+M y_c(t)+N d(t) \]

更普遍的,将扰动模型写作:

\[ \begin{aligned}& \dot{x}_s(t)=A_s x_s(t) \\& d(t)=C_s x_s(t)\end{aligned} \]

由此可以构建增广模型:

\[ \begin{gathered}\binom{\dot{x}(t)}{\dot{x}_s(t)}=\left(\begin{array}{cc}A & B^{\prime} C_s \\0 & A_s\end{array}\right)\binom{x(t)}{x_s(t)}+\binom{B}{0} u(t)=A_a x_a(t)+B_a u(t) \\y(t)=\left(\begin{array}{ll}C & D^{\prime} C_s\end{array}\right)\binom{x(t)}{x_s(t)}=C_a x_a(t)\end{gathered} \]

得到如下形式的控制律:

\[ u(t)=-\left(\begin{array}{ll} K & K_s \end{array}\right)\binom{x(t)}{x_s(t)}+e(t)=-K_a x_a(t)+e(t) \]


不可测量的情况

在这种情况下,观测器不仅要重构状态 \(x(t)\) ,还要重构附加的状态变量 \(x_s(t)\) 。构造的观测器闭环系统为:

\[ \begin{aligned}& \binom{\dot{\hat{x}}(t)}{\dot{\hat{x}}_s(t)}=\left(\begin{array}{cc}A & B^{\prime} C_s \\0 & A_s\end{array}\right)\binom{\hat{x}(t)}{\hat{x}_s(t)}+\binom{B}{0} u(t)+\binom{L}{L_s}\left(y(t)-\left(\begin{array}{ll}C & D^{\prime} C_s\end{array}\right)\binom{\hat{x}(t)}{\hat{x}_s(t)}\right) \\& \quad=A_a \hat{x}_a(t)+B_a u(t)+L_a\left(y(t)-C_a \hat{x}_a(t)\right)\end{aligned} \]

构建的控制律:

\[ u(t)=-K \hat{x}(t)+M y_c(t)+N C_s\hat{x}_s(t) = -\left(\begin{array}{ll} K & K_s \end{array}\right)\binom{\hat{x}(t)}{\hat{x}_s(t)}+e(t)=-K_a \hat{x}_a(t)+e(t) \]

CE4.3 电磁悬挂系统


讨论电磁悬挂系统在状态不可测量的其概况。

根据之前章节的讨论,系统状态方程为:

\[ \begin{aligned} \frac{d}{d t}\left(\begin{array}{c} i_1(t) \\ z(t) \\ \dot{z}(t) \end{array}\right) & =\left(\begin{array}{ccc} -1000 & 0 & 0 \\ 0 & 0 & 1 \\ 20 & 1600 & 0 \end{array}\right)\left(\begin{array}{c} \dot{i}_1(t) \\ z(t) \\ \dot{z}(t) \end{array}\right)+\left(\begin{array}{c} 100 \\ 0 \\ 0 \end{array}\right) u_1(t) \\ V_z(t) & =\left(\begin{array}{lll} 0 & 4000 & 0 \end{array}\right)\left(\begin{array}{l} \dot{i}_1(t) \\ z(t) \\ \dot{z}(t) \end{array}\right) \end{aligned} \]

系统完全可观测。观测器模型为:

\[ \begin{aligned}\frac{d}{d t}\left(\begin{array}{l} \hat{i}_1(t) \\ \hat{z}(t) \\ \hat{z}(t) \end{array}\right)=&\left(\begin{array}{ccc} -1000 & 0 & 0 \\ 0 & 0 & 1 \\ 20 & 1600 & 0 \end{array}\right)\left(\begin{array}{c} \hat{1}_1(t) \\ \hat{z}(t) \\ \hat{\dot{z}}(t) \end{array}\right)+\left(\begin{array}{c} 100 \\ 0 \\ 0 \end{array}\right) u_1(t)+\left(\begin{array}{l} l_1 \\ l_2 \\ l_3 \end{array}\right)\left(V_z(t)-\left(\begin{array}{lll} 0 & 4000 & 0 \end{array}\right)\left(\begin{array}{l} \hat{1}_1(t) \\ \hat{z}(t) \\ \hat{\dot{z}}(t) \end{array}\right)\right)\\=& \left(\begin{array}{ccc} -1000 & -4000 l_1 & 0 \\ 0 & -4000 l_2 & 1 \\ 20 & 1600-4000 l_3 & 0 \end{array}\right)\left(\begin{array}{c} \hat{i}_1(t) \\ \hat{z}(t) \\ \hat{z}(t) \end{array}\right)+\left(\begin{array}{c} 100 \\ 0 \\ 0 \end{array}\right) u_1(t)+\left(\begin{array}{l} l_1 \\ l_2 \\ l_3 \end{array}\right) V_z(t) \end{aligned} \]

在之前的选择中,我们使用了\(\left\{\begin{array}{l}-1000\\\xi=0,6 \\\omega_0=130\end{array}\right.\)的参数。我们希望观察器能有3倍的速度。回忆收敛速度与特征值的负实部的绝对值呈正比。而复共轭特征值通常表示为\(\lambda_{1,2}=-\zeta \omega_0 \pm j \omega_0 \sqrt{1-\zeta^2}\),参考一年前的笔记如下:

因此,需选择\(\left\{\begin{array}{l}-3000\\\xi=0,6 \\\omega_0=390\end{array}\right.\)对应的特征值:

\[ |\lambda \mathbf{I}-(\mathbf{A}-\mathbf{L} \mathbf{C})|=(\lambda+3000)\left(\lambda^2+2 \times 0,6 \times 390 \lambda+390^2\right)\Rightarrow \begin{aligned} & l_1=17100 \\ & l_2=0,617 \\ & l_3=-227,6 \end{aligned} \]

CE3.4 海洋种群


海洋种群

🪺

我们考虑以下数学模型:

\[ \left\{\begin{array}{l}\dot{z}_1(t)=-z_1(t)+z_1(t) z_2(t) \\\dot{z}_2(t)=-z_2(t)-z_1(t) z_2(t)+z_2(t) z_3(t) \\\dot{z}_3(t)=-z_3(t)-z_2(t) z_3(t)+v(t)\end{array}\right. \]

其中, \(z_1\) 代表主要捕食者的种群(例如虎鲸), \(z_2\) 代表中间捕食者/猎物种群(例如海豹), \(z_3\) 代表最终猎物种群(例如草食性鱼类), \(v\) 代表这些猎物的食物(例如藻类)。

考虑\(v(t)=\bar{v}\)时的平衡点:

\[ \left\{\begin{array}{l}0=-\bar{z}_1+\bar{z}_1 \bar{z}_2=\bar{z}_1\left(\bar{z}_2-1\right) \\0=-\bar{z}_2-\bar{z}_1 \bar{z}_2+\bar{z}_2 \bar{z}_3=\bar{z}_2\left(-1-\bar{z}_1+\bar{z}_3\right) \\0=-\bar{z}_3-\bar{z}_2 \bar{z}_3+\bar{v}=\bar{z}_3\left(1+\bar{z}_2\right)+\bar{v}\end{array}\right. \]

特别考虑以下爱平衡点:

\[ \left(\bar{z}_1, \bar{z}_2, \bar{z}_3\right)=\left(\frac{\bar{v}-2}{2}, 1, \frac{\bar{v}}{2}\right) \]

在平衡点附近线性化模型:

\[ \dot{x}(t)=\left(\begin{array}{ccc}-1+\bar{z}_2 & \bar{z}_1 & 0 \\-\bar{z}_2 & -1-\bar{z}_1+\bar{z}_3 & \bar{z}_2 \\0 & -\bar{z}_3 & -1-\bar{z}_2\end{array}\right) x(t)+\left(\begin{array}{l}0 \\0 \\1\end{array}\right) u(t) \]

其中\(x=z-\bar{z}, u=v-\bar{v}\),设\(\bar{v}=3\),有:

\[ \dot{x}(t)=A x(t)+B u(t)=\left(\begin{array}{ccc}0 & 0.5 & 0 \\-1 & 0 & 1 \\0 & -1.5 & -2\end{array}\right) x(t)+\left(\begin{array}{l}0 \\0 \\1\end{array}\right) u(t) \]

\(spectre(A) = [-1,-0.5 \pm 0.866 i]\),因此系统是稳定的。同时,控制矩阵:

\[ rang(C)=rang(\left[\begin{array}{lll} B & A B & A^2 B \end{array}\right])=rang(\left(\begin{array}{ccc} 0 & 0 & 0.5 \\ 0 & 1 & -2 \\ 1 & -2 & 2.5 \end{array}\right)) = 3 \]

因此,系统可控。


首先加色所有变量都可测量,通过状态反馈控制,系统闭环模型为:

\[ \dot{x}(t)=(A-B K) x(t)+B M e(t)=\left(\begin{array}{ccc}0 & 0.5 & 0 \\-1 & 0 & 1 \\-k_1 & -1.5-k_2 & -2-k_3\end{array}\right) x(t)+\left(\begin{array}{c}0 \\0 \\M\end{array}\right) y_c(t) \]

其特征多项式:\(\operatorname{det}(\lambda I-(A-B K))=\lambda^3+\left(2+k_3\right) \lambda^2+\left(2+k_2\right) \lambda+\left(1+0.5 k_1+0.5 k_3\right)\)

格定义两个复共轭特征值,其阻尼比为 \(\xi=0.46\) ,固有频率为 \(\omega_0=0.354\) 。第三个特征值选择在前两个特征值的左侧较远处,以避免改变所选的瞬态行为,优先选择 \(-2\) 。其特征多项式:

\[ \operatorname{det}(\lambda I-(A-B K))=\left(\lambda^2+2 \xi \omega_0 \lambda+\omega_0^2\right)(\lambda+2) \]

得到:

\[ \left\{\begin{array}{l}k_1=4 \omega_0^2-2 \xi \omega_0-2 \simeq-1,825 \\k_2=\omega_0^2+4 \xi \omega_0-2 \simeq-1,224 \\k_3=2 \xi \omega_0 \simeq 0,326 .\end{array}\right. \]

为了确定增益 \(M\) ,我们考察了应用于系统输入的幅值为 \(y_{c 0}\) 的阶跃信号的影响。在这种输入下,线性化系统的平衡点 \(\bar{x}\) 由以下方程给出:

\[ \left(\begin{array}{l}0 \\0 \\0\end{array}\right)=\left(\begin{array}{ccc}0 & 0.5 & 0 \\-1 & 0 & 1 \\-k_1 & -1.5-k_2 & -2-k_3\end{array}\right)\left(\begin{array}{l}\bar{x}_1 \\\bar{x}_2 \\\bar{x}_3\end{array}\right)+\left(\begin{array}{c}0 \\0 \\M\end{array}\right) y_{c 0} \]

即:

\[ \begin{gathered}\bar{x}_2=0 \\\bar{x}_1=\bar{x}_3=\frac{M y_{c 0}}{2+k_1+k_3}\end{gathered} \]

我们可以通过选择\(M=2+k_1+k_3 \approx 0.501\)


现在假设只有草食性鱼类种群 \(x_3\)是可测量的,重建模型为:

\[ \dot{\hat{x}}(t)=A \hat{x}(t)+B u(t)+L(y(t)-C \hat{x}(t)) \]

假设只有\(x_3\)是可观察的:

\[ y(t)=C x(t)=\left(\begin{array}{lll}0 & 0 & 1\end{array}\right) x(t) \]

验证其可观察性:

\[ \mathcal{O}=\left[\begin{array}{c} C \\ C A \\ C A^2 \end{array}\right]=\left(\begin{array}{ccc} 0 & 0 & 1 \\ 0 & -1,5 & -2 \\ 1,5 & 3 & 2,5 \end{array}\right) \]

满秩。因此新的模型写作:

\[ \begin{aligned} \dot{\hat{x}}(t) & =A \hat{x}(t)+B u(t)+L(y(t)-C \hat{x}(t)) \\ & =\left(\begin{array}{ccc} 0 & 0,5 & 0 \\ -1 & 0 & 1 \\ 0 & -1,5 & -2 \end{array}\right) \hat{x}(t)+\left(\begin{array}{l} 0 \\ 0 \\ 1 \end{array}\right) u(t)+\left(\begin{array}{c} \ell_1 \\ \ell_2 \\ \ell_3 \end{array}\right)\left(y(t)-\hat{x}_3(t)\right) \end{aligned} \]

其中增益向量 \(L=\left(\ell_1, \ell_2, \ell_3\right)^{\top}\) 需要确定。我们为观测器选择两个复共轭特征值,阻尼比为 \(\xi^{\prime}=0.7\) (以限制误差估计中的瞬态振荡),固有频率为原来计算中选取频率的 \(2.5\) 倍(以确保估计比控制更快地进行),即 \(\omega_1=0.885\) ,并在 -2 处再添加一个特征值(为什么不选择\(2.5\)倍呢)。计算得到:

\[ \operatorname{det}(\lambda I-(A-L C))=\lambda^3+\left(2+\ell_3\right) \lambda^2+\left(2-1.5 \ell_2\right) \lambda+\left(1+1.5 \ell_1+0.5 \ell_3\right) = \operatorname{det}(\lambda I-(A-L C))=\left(\lambda^2+2 \xi^{\prime} \omega_1 \lambda+\omega_1^2\right)(\lambda+2) \]

得到解:

\[ \left\{\begin{array}{l}\ell_1=\frac{2 \omega_1^2-\xi^{\prime} \omega_1-1}{1.5} \approx-0.0354 \\\ell_2=\frac{-\omega_1^2-4 \xi^{\prime} \omega_1+2}{1.5} \approx-0.8408 \\\ell_3=2 \xi^{\prime} \omega_1 \approx 1.239\end{array}\right. \]

估计控制律为:\(u(t)=-K \hat{x}(t)+M y_c(t) ==-\left(\begin{array}{lll} -1.825 & -1.224 & 0.326) \hat{x}(t)+0.501 y_c(t) \end{array}\right.\)


上图展示了模拟的结果。在过渡阶段,估计状态变量和实际状态变量之间存在显著差异。然而,由于估计误差的收玫(即矩阵 \(A-L C\) 的稳定性),这种差异很快就消失了。注意,如果我们将观测器初始化为 \(\hat{x}(0)=x(0)\) ,那么所获得的响应将与测量结果完全一致(实线图),因为正如我们所见,估计误差将为零。


考虑非线性系统,重新控制律:

\[ \begin{gathered} v(t)=-K \hat{x}(t)+M y(t)+\bar{v} \\ \hat{x}(t)=A \hat{x}(t)+B(v(t)-\bar{v})+L\left(z_3(t)-\bar{z}_3-C \hat{x}(t)\right) \end{gathered} \]

为了保持线性化有效性的邻域内,我们考虑幅值较小的设定值阶跃信号。可以观察到,使用状态反馈控制 (虚线曲线) 的性能与线性化系统的预测几乎相同。这是因为较小的设定值阶跃信号使得状态 \(z\) 保持在平衡点 \(\bar{z}\) 附近,足够接近于线性化系统的行为。

C4.2 优化方法 Approche optimale


对于有噪音的情况:

\[ \begin{aligned} & \dot{x}(t)=A x(t)+B u(t)+v(t) \\ & y(t)=C x(t)+w(t) . \end{aligned} \]

其中\(v(t), w(t)\) 是均值为零的高斯白噪声且不相关:

\[ \left\{\begin{array}{l} \mathbb{E}\{v(t)\}=0 \\ \mathbb{E}\left\{v(t) v(\tau)^{\top}\right\}=V \delta(t-\tau) \end{array}, \quad\left\{\begin{array}{l} \mathbb{E}\{w(t)\}=0 \\ \mathbb{E}\left\{w(t) w(\tau)^{\top}\right\}=W \delta(t-\tau), \end{array}\right.\right. \]

此时观察器的结构为:

\[ \dot{\hat{x}}(t)=A \hat{x}(t)+B u(t)+L(y(t)-C \hat{x}(t)), \]

其误差\(\tilde{x}=x-\hat{x}\)和误差动态方程:

\[ \begin{aligned} \dot{\tilde{x}}(t) & =\dot{x}(t)-\dot{\hat{x}}(t) \\ & =A x(t)+B u(t)+v(t)-(A \hat{x}(t)+B u(t)+L(y(t)-C \hat{x}(t))) \\ & =A \tilde{x}(t)+v(t)-L(C x(t)+w(t)-C \hat{x}(t)),\\&=(A-L C) \tilde{x}(t)+v(t)-L w(t) \end{aligned} \]

为此,如何在存在噪声和模型不确定性的情况下确定增益 \(L\),以实现最优且无偏估计?


目标:从 \(u(t)\)\(y(t)\) 重构 \(x(t)\)

  • 目标\(1\):无偏估计estimation sans biais\(\dot{\tilde{x}}(t)=\dot{x}(t)-\dot{\hat{x}}(t) \text { t.q. } E\{\boldsymbol{\varepsilon}(t)\}=0\)
  • 目标\(2\):最小化估计误差的方差,\(\boldsymbol{\Sigma}(t)=E\left\{\tilde{x}(t) \tilde{x}(t)^T\right\}\)

考虑目标1:

\[ \begin{aligned} \mathbb{E}\{\tilde{x}(t)\} & =e^{(A-L C) t} \mathbb{E}\{\tilde{x}(0)\}+\int_0^t e^{(A-L C)(t-\tau)} \mathbb{E}\{v(\tau)-L w(\tau)\} d \tau \\ & =e^{(A-L C) t} \mathbb{E}\{\tilde{x}(0)\} . \end{aligned} \]

因此,有两种可行方式:

\[ \begin{aligned}& E\{\boldsymbol{\varepsilon}(t)\}=0 \ \forall t, \text { si } E\{\boldsymbol{\varepsilon}(0)\}=0 \\& E\{\boldsymbol{\varepsilon}(t)\} \underset{t \rightarrow \infty}{\rightarrow} 0, \text { si } A-LC \ \text{stable }\end{aligned} \]

前者要求\(\hat{x}(0)=E\{x(0)\}\),后者要求\(A-LC\)稳定。


考虑目标2:

为了满足目标2,应该选择\(L\),使其最小化估计误差的方差:

\[ \Sigma(t)=\mathbb{E}\left\{\tilde{x}(t) \tilde{x}(t)^{\top}\right\} \]

C4.2.1 卡尔曼滤波器 Filtre de Kalman


定理

允许最小化估计误差方差 \(\Sigma(t)=\mathbb{E}\left\{\tilde{x}(t) \tilde{x}(t)^{\top}\right\}\) 的矩阵 \(L(t)\) 在任何时刻 \(t \geq 0\) 给出为

\[ L(t)=\Sigma(t) C^{\top} W^{-1} \]

其中 \(\Sigma(t)\) 由Riccati微分方程的解给出

\[ \dot{\Sigma}(t)=A \Sigma(t)+\Sigma(t) A^{\top}+V-\Sigma(t) C^{\top} W^{-1} C \Sigma(t) \]

初始条件为 \(\Sigma(0)=\mathbb{E}\left\{\tilde{x}(0) \tilde{x}(0)^{\top}\right\}\)


有两点需要注意:

  • 尽管系统是时间不变的,观测器使用的矩阵 \(L\) 不是常数。它依赖于估计误差的方差 \(\Sigma\) ,该方差随时间根据动态方程演变。
  • 在实际操作中,无法知道 \(\mathbb{E}\left\{\tilde{x}(0) \tilde{x}(0)^{\top}\right\}\) (这需要知道初始的观测误差,因此也需要知道初始状态本身)。因此,不能按照定理中的方式初始化 \(\Sigma\) 。初始值 \(\Sigma(0)\) 被作为一个参数来调整,以改善得到的性能。

Riccati 微分方程的解法

将Riccati方程写作:

\[ \left[\begin{array}{c}\dot{X}(t) \\\dot{Y}(t)\end{array}\right]=\left[\begin{array}{cc}-A^{\top} & C^{\top} W^{-1} C \\V & A\end{array}\right]\left[\begin{array}{l}X(t) \\Y(t)\end{array}\right] \]

其中,\(\Sigma(t)=Y(t) X(t)^{-1}\),以及初始条件:\(X(0)=I ,Y(0)=\mathbb{E}\left\{\tilde{x}(0) \tilde{x}(0)^{\top}\right\}\)

因为:

\[ \dot{\Sigma}(t)=\dot{Y}(t) X(t)^{-1}-Y(t) X(t)^{-1} \dot{X}(t) X(t)^{-1} \]

带入Riccati方程得到:

\[ \begin{aligned}\dot{\Sigma}(t) & =(V X(t)+A Y(t)) X(t)^{-1}-Y(t) X(t)^{-1}\left(-A^{\top} X(t)+C^{\top} W^{-1} C Y(t)\right) X(t)^{-1} \\& =V+A Y(t) X(t)^{-1}+Y(t) X(t)^{-1} A^{\top}-Y(t) X(t)^{-1} C^{\top} W^{-1} C Y(t) X(t)^{-1} \\& =V+A \Sigma(t)+\Sigma(t) A^{\top}-\Sigma(t) C^{\top} W^{-1} C \Sigma(t),\end{aligned} \]

从而证明上述写法合理。求解得:

\[ \left[\begin{array}{l} X(t) \\ Y(t) \end{array}\right]=e^{H t}\left[\begin{array}{c} X(0) \\ Y(0) \end{array}\right]=e^{H t}\left[\begin{array}{c} I \\ \mathbb{E}\left\{\tilde{x}(0) \tilde{x}(0)^{\top}\right\} \end{array}\right]=e^{H t}\left[\begin{array}{c} I \\ \Sigma(0) \end{array}\right] \]

其中,

\[ H= \left[\begin{array}{cc} -A^{\top} & C^{\top} W^{-1} C \\ V & A \end{array}\right] \]

CE3.5 简单的示例


设:

\[ \begin{gathered}\left\{\begin{array}{l}\dot{x}(t)=u+V \\y=x+W\end{array}\right. ,\quad V=W=1\end{gathered} \]

有:

\[ H= \left[\begin{array}{cc} -A^{\top} & C^{\top} W^{-1} C \\ V & A \end{array}\right] = \left[\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array}\right] \]

又有:

\[ \exp \left(\left[\begin{array}{ll}0 & 1 \\1 & 0\end{array}\right] t\right)=\left[\begin{array}{ll}\cosh t & \sinh t \\\sinh t & \cosh t\end{array}\right] \]

最终得到:

\[ \Sigma(t)=\frac{1-e^{-2 t}+\left(1+e^{-2 t}\right) \Sigma(0)}{1+e^{-2 t}+\left(1-e^{-2 t}\right) \Sigma(0)} \]


Propriété de l'innovation

对于上述系统,假设\(\dot{\hat{x}}(t)=A \hat{x}(t)+B u(t)+L(y(t)-C \hat{x}(t))\),称:

\[ \tilde{y}(t)=y(t)-C \hat{x}(t) \]

为创新(innovation),是一个白噪音。

C4.2.2 Kalman滤波器的渐进行为


卡尔曼滤波器中使用的矩阵随时间变化,即使所考虑的系统是时间不变的。我们研究\(L(t)\)收敛到一个接近于常数值得条件,并选择该常数增益。

假设\(V=V_0 V_0^{\top}\) ,那么,使用渐近增益\(L\) 的卡尔曼滤波器稳定的一个充分条件是对偶对\((A, V_0)\)是可稳定的且对偶对\((A, C)\)是可检测的。


如果满足上述条件,则有\(\Sigma(t)\)收敛到下述方程得解\(\Sigma\)

\[ A \Sigma+\Sigma A^{\top}+V-\Sigma C^{\top} W^{-1} C \Sigma=0 \]

若取\(L=\Sigma C^{\top} W^{-1}\),则常生的估计误差将最终收敛到\(0\)

最终,如果\((A, V_0)\)是可控的,则\(\Sigma\)是正定的。这一渐近值并不保证最小化估计误差的方差。

C4.2.3 与LQ控制的对偶性


回忆LQ控制:

\[ K=R^{-1} B^{\top} P,A^{\top} P+P A-P B R^{-1} B^{\top} P+Q=0 \]

可以发现,LQ控制与Kalman滤波器对偶。

\[ L=\Sigma C^T W^{-1},A \Sigma+\Sigma A^{\top}+V-\Sigma C^{\top} W^{-1} C \Sigma=0 \]

有:

\[ \mathbf{K} \leftrightarrow \mathbf{L}^T, \mathbf{A} \leftrightarrow \mathbf{A}^T, \mathbf{B} \leftrightarrow \mathbf{C}^T, \mathbf{C} \leftrightarrow \mathbf{B}^T, \mathbf{P} \leftrightarrow \mathbf{\Sigma}, \mathbf{Q} \leftrightarrow \mathbf{V}, \mathbf{R} \leftrightarrow \mathbf{W} \]

C4.2.4 Kalman滤波器的调节参数


  • \(V, W\) 是卡尔曼滤波器的调节参数。
    • \(V\) :对系统模型的信任度(状态方程)
    • \(W\) :对测量的信任度
  • 实际上,关煡在于 \(V\)\(W\) 之间的"比例"。
    • 如果 \(V\) 相对于 \(W\) 较小:更信任模型
      • \(L\) 较小
      • 估计结果较平滑
      • 估计速度较慢
    • 如果 \(V\) 相对于 \(W\) 较大:更信任测量
      • \(L\) 较大
      • 估计结果较噪声化
      • 估计速度较快

C4.3 高斯线性二次型控制 (LQG)


考虑使用卡尔曼滤波器来控制只能部分状态可测的系统。

\[ \begin{gathered} \dot{x}(t)=A x(t)+B u(t)+v(t) \\ y(t)=C x(t)+w(t), \end{gathered} \]

其中\(y\)可测。由于存在噪音,状态不会收敛到0,为了避免可能的\(\int_0^{\infty}\left(x(t)^{\top} Q x(t)+u(t)^{\top} R u(t)\right) d t\),使用以下准则:

\[ J=\lim _{T \rightarrow+\infty} \mathbb{E}\left\{\frac{1}{T} \int_0^T\left(x(t)^{\top} Q x(t)+u(t)^{\top} R u(t)\right) d t\right\} \]

相比于之前讨论的情况,x不再可测。唯一的观测器是:

\[ \hat{x}(t)=A \hat{x}(t)+B u(t)+L(t)(y(t)-C \hat{x}(t)) \]

\(Q=Q_0 Q_0^{\top}\),假定\((A,B)\)可稳定,\((A,Q_0)\)可检测。则最小化二次型准则的线性控制位:

\[ \begin{aligned} & u(t)=-K \hat{x}(t) \\ & K=R^{-1} B^{\top} P, \end{aligned} \]

其中\(P\)是Riccati方程的半正定对称解:

\[ A^{\top} P+P A-P B R^{-1} B^{\top} P=-Q \]

\(\hat x\)由观测器给出,且\(L(t)=\Sigma(t) C^{\top} W^{-1}\),满足:

\[ \dot{\Sigma}(t)=\Sigma(t) A^{\top}+A \Sigma(t)+V(t)-\Sigma(t) C^{\top} W^{-1} C \Sigma(t) \]

在这种情况下,我们将先确定\(\Sigma\)\(L\),然后再确认\(P\)\(K\)

C2.6 飞机的水平运动问题


我们曾经讨论过这个例子,因此不再赘述。我们假设:

  • 侧滑角 \(\beta\) 无法测量
  • 横滚率 \(p\) 、偏航率 \(r\) 和横滚角 \(\phi\) 的测量值受到噪声影响,分别记为 \(w_1 、 w_2\)\(w_3\)
  • 系统受到附加在控制输入 \(\delta_a\)\(\delta_b\) 上的噪声 \(b_1\)\(b_2\) 的激励。

在线性化后,系统为:

\[ \begin{gathered}\dot{x}(t)=A x(t)+B(u(t)+b(t)) \\y(t)=C x(t)+w(t)\end{gathered} \]

其中:

\[ x=(\beta, p, r, \phi)^{\top}, \quad u=\left(\delta_a, \delta_b\right)^{\top}, \quad b=\left(b_1, b_2\right)^{\top}, \quad w=\left(w_1, w_2, w_3\right)^{\top}, \]

\[ A=\left(\begin{array}{cccc}-0,746 & 0,006 & -1 & 0,0369 \\-12,9 & -0,746 & 0,387 & 0 \\4,31 & 0,024 & -0,174 & 0 \\0 & 1 & 0 & 0\end{array}\right), B=\left(\begin{array}{ccc}0,0012 & 0,0092 \\6,05 & 0,952 \\-0,416 & -1,76 \\0 & 0\end{array}\right), \quad C=\left(\begin{array}{llll}0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1\end{array}\right) . \]

对于噪音,我们假设她们彼此不相关:

\[ \begin{aligned} \mathbb{E}\left\{b(t) b(t)^{\top}\right\} & =\left(\begin{array}{cc} \alpha_1 & 0 \\ 0 & \alpha_2 \end{array}\right) \\ \mathbb{E}\left\{w(t) w(t)^{\top}\right\} & =W=\left(\begin{array}{ccc} \gamma_1 & 0 & 0 \\ 0 & \gamma_2 & 0 \\ 0 & 0 & \gamma_3 \end{array}\right) . \end{aligned} \]

假设\(v(t)=B b(t)\),噪音\(v\)的协方差矩阵位:

\[ V=\mathbb{E}\left\{v(t) v(t)^{\top}\right\}=\mathbb{E}\left\{B b(t) b(t)^{\top} B^{\top}\right\}=B\left(\begin{array}{cc} \alpha_1 & 0 \\ 0 & \alpha_2 \end{array}\right) B^{\top} \]

故有观察器:

\[ \dot{\hat{x}}(t)=A \hat{x}(t)+B u(t)+L(y(t)-C \hat{x}(t)) \]

其中:

\[ L=\Sigma C^{\top} W^{-1} \]

\(\Sigma\)是Riccati方程在稳态下的解:

\[ A \Sigma+\Sigma A^{\top}+V-\Sigma C^{\top} W^{-1} C \Sigma=0 \]

在验证\(\left(A, V_0\right)\),其中\(V=V_0 V_0^{\top}, \quad V_0=B\left(\begin{array}{cc} \sqrt{\alpha_1} & 0 \\ 0 & \sqrt{\alpha_2} \end{array}\right)\),是可控的之后,可以对方程进行数值求解,假设\(\gamma_1=\gamma_2=10^{-3} 、 \gamma_3=10^{-4}\),假设\(\gamma_1=\gamma_2=10^{-3} 、 \gamma_3=10^{-4}\),得到:

\[ \Sigma=\left(\begin{array}{cccc} 4.9352 & -2.8245 & -5.7488 & -0.0512 \\ -2.8245 & 185.7161 & -13.5804 & 8.2863 \\ -5.7488 & -13.5804 & 49.1163 & -0.2865 \\ -0.0512 & 8.2863 & -0.2865 & 3.1142 \end{array}\right) \times 10^{-4} \]

遂有:

\[ L=\Sigma C^{\top} W^{-1}=\left(\begin{array}{ccc} -0.2824 & -0.5749 & -0.0512 \\ 18.5716 & -1.3580 & 8.2863 \\ -1.3580 & 4.9116 & -0.2865 \\ 0.8286 & -0.0286 & 3.1142 \end{array}\right) \]

验证其正定性,得:

\[ \operatorname{spectre}(A-L C)=\{-19,1630,-1,3149,-4,5727,-3,2128\} \]

符合。

进而将进行LQ控制,假设:

\[ Q=\left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 10 \end{array}\right) \quad \text { et } \quad R=\left(\begin{array}{ll} 3 & 0 \\ 0 & 1 \end{array}\right) \]

有:

\[ Q_0=\left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \sqrt{10} \end{array}\right) \]

可验证\((A,Q_0)\)可观测,继而考虑:

\[ \begin{aligned} u(t) & =-K \hat{x}(t) \\ K & =R^{-1} B^{\top} P \\ \dot{\hat{x}}(t) & =A \hat{x}(t)+B u(t)+L(y(t)-C \hat{x}(t)) \end{aligned} \]

得:

\[ A^{\top} P+P A-P B R^{-1} B^{\top} P=-Q\\P=\left(\begin{array}{cccc} 4,0543 & -0,9381 & -1,0488 & -1,7819 \\ -0,9381 & 0,3440 & 0,2830 & 0,9135 \\ -1,0488 & 0,2830 & 0,8213 & 0,2837 \\ -1,7819 & 0,9135 & 0,2837 & 4,1934 \end{array}\right)\\K=R^{-1} B^{\top} P=\left(\begin{array}{cccc} -1,745 & 0,654 & 0,456 & 1.802 \\ 0,99 & -0,179 & -1,186 & 0.354 \end{array}\right) \]