CN104503428A - 一种民机飞控系统抗干扰时变故障诊断方法 - Google Patents

Abstract

一种民机飞控系统抗干扰时变故障诊断方法。其包括按顺序进行的下列步骤：建立包含外部干扰和内部故障向量的民机飞控系统非线性数学模型；设定故障估计算法的约束条件；设计状态观测器；设计自适应故障估计算法。本发明优点：①不是基于理想解析模型的，在设计阶段，外部干扰已被事先考虑进来，因此提高了后续故障诊断准确率。②可同时获得系统状态变量和故障的具体形态。③对故障类型不敏感。④外部干扰和系统内部组件故障均可以是快速时变的(非常值或慢时变类型)。⑤本方法是基于非线性系统模型的，而不是传统的基于线性模型的故障诊断方法。⑥算法规模小、运算简单，易于工程实现。本技术可进一步推广应用于一般非线性控制系统的故障诊断。

Description

一种民机飞控系统抗干扰时变故障诊断方法

技术领域

本发明属于故障诊断与容错控制技术领域，特别是涉及一种民机飞控系统抗干扰时变故障诊断方法。

背景技术

近年来，受民航业发展驱动，民用飞机自动飞行控制系统(下文简称民机飞控系统)的复杂度、性能和自动化程度持续提高，这就要求民机飞控系统较之以前必须具备更好的可靠性和安全性。目前，民机飞控系统可靠性问题已经成为飞机设计、制造领域中最关键、最亟待解决的技术问题之一。增加系统可靠性最直接的一种方法是提高系统组成部件的质量，如执行器、控制器或传感器，但是，这样会大幅度增加设备的支出成本，代价昂贵，而且即使这样也很难避免绝对不会发生故障。基于此，民机飞控系统的故障诊断问题变得愈发重要，目前大量科研单位的技术工作人员正致力于该项研究。

民机飞控系统故障诊断方法中具有代表性的是Alcorta-Garcia Efrain等在2011年提出的小幅度、振荡型故障的早期探测方法。该方法通过对A380机型飞行控制系统组成的深入分析，建立了理想的解析非线性模型，然后设计出故障诊断算法，通过缩小探测阈值，实现了小幅值故障的早期探测并给出故障告警，还通过数据仿真实验证明了所提出算法的有效性。但同时也存在几个问题，如：(1)所用模型是理想的解析数学模型，没有考虑外部干扰的影响，很容易造成最终诊断结果中干扰与故障的混淆，即外部干扰并不是故障，因此此种诊断方法不能将干扰与故障很好地区分开来。(2)诊断结果的输出形式为超限警告，即输出故障提醒信息，但无法给出所发生故障的具体形式。(3)所提出的故障诊断算法仅仅适用于“振荡型”故障，对于其他类型的故障，如“斜坡型”、“抛物线型”等类型的故障无效。(4)所提出的故障诊断算法复杂、规模和计算量大，工程实现时的难度较大。

此外，民机飞控系统的故障诊断一般只考虑对故障向量的估计与检测，但是在发生故障时系统状态变量的运动状态对后续的故障原因分析以及容错控制也是非常重要的。L.Lavigne将一种改良的卡尔曼滤波器用于故障和系统状态变量的同时估计，但是这种方法仅适用于常值或慢时变类型故障的探测，对于快速时变故障的探测性能极剧下降，因此容易出现漏报或虚警现象。

发明内容

为了解决上述问题，本发明的目的在于提供一种能够在给出故障具体形态的同时有效地将系统状态变量也估计出来，并且算法规模小、运算简单、易于工程实现的民机飞控系统抗干扰时变故障诊断方法。

为了达到上述目的，本发明提供的民机飞控系统抗干扰时变故障诊断方法包括按顺序进行的下列步骤：

(1)建立同时包含外部干扰和内部故障向量的民机飞控系统非线性数学模型；

(2)设定故障估计算法的约束条件；

(3)设计状态观测器；

(4)设计自适应故障估计算法。

在步骤(1)中，所述的民机飞控系统非线性数学模型是：

$\left\{\begin{array}{c}\stackrel{.}{x}\left(t\right)=\mathrm{\φ}(x,t)+B[u\left(t\right)+\mathrm{\γ}\left(t\right)]+\mathrm{Ed}\left(t\right)\\ y\left(t\right)=\mathrm{Cx}\left(t\right)\end{array}\right.---\left(1\right)$

其中，表示状态向量；表示控制输入；表示输出向量；表示系统受到的外部时变干扰；φ(x,t)是一个关于x(t)的非线性函数；表示未知的待诊断时变故障；B，C和E均为适当维的常值矩阵。

在步骤(2)中，所述的故障估计算法的约束条件是：

①对于状态向量x(t)，非线性函数φ(x,t)满足Lipschitz条件，即存在一个常数δ₁>0使得如下不等式成立：

||φ(x₁)-φ(x₂)||≤δ₁||x₁-x₂||   (2)

x(t)缩写成x；

②外部时变干扰d的能量是有界的，即：

||Ed||≤η₀   (3)

其中，η₀≥0表示系统能够承受的最大干扰量；

③未知的待诊断的时变故障向量γ的能量和变化频率均有界，即：

||γ||≤η₁

                                               (4) 

$\left|\right|\stackrel{.}{\mathrm{\γ}}\left|\right|\≤{\mathrm{\η}}_2$

其中，η₁≥0和η₂≥0分别表示发生故障的最大幅值和最快频率；

④状态向量x的能量和变化频率均有界，即：

||x||≤η₃

                                               (5) 

$\left|\right|\stackrel{.}{x}\left|\right|\≤{\mathrm{\η}}_4$

其中，η₃≥0和η₄≥0分别表示状态变量的最大幅值和最大变化速率。

在步骤(3)中，所述的状态观测器的设计方法是：

利用步骤(1)得到的民机飞控系统非线性数学模型和步骤(2)得到的四个约束条件，设计如下式(6)所示的状态观测器，用于完成外部干扰和内部时变故障同时存在时系统状态变量x的估计；

$\left\{\begin{array}{c}\stackrel{\stackrel{.}{^}}{x}=\mathrm{\φ}\left(\hat{x}\right)+B[u+\widehat{\mathrm{\γ}}]+L(y-\hat{y})\\ \hat{y}=C\hat{x}\end{array}\right.---\left(6\right)$

其中，表示对系统状态变量x的估计值；表示对时变故障γ的估计值； 表示对系统输出y的估计值；L是一个可调节增益矩阵。

在步骤(4)中，所述的自适应故障估计算法的设计方法是：

结合步骤(1)得到的民机飞控系统非线性数学模型和步骤(3)得到的状态观测器，设计如下式(7)所示的自适应故障估计算法，用于完成对系统内部时变故障γ的精确估计：

$\stackrel{\stackrel{.}{^}}{\mathrm{\γ}}={\mathrm{\ΓB}}^{T}P\tilde{x}+M\stackrel{\stackrel{.}{~}}{x}---\left(7\right)$

其中，表示状态观测器的观测误差；表示误差估计的误差；  $\mathrm{\Γ}=\mathrm{diag}\{{\mathrm{\Γ}}_1,{\mathrm{\Γ}}_2,...,{\mathrm{\Γ}}_r\},\∀{\mathrm{\Γ}}_i>0,i\∈[1,r]$ 和 $M=\mathrm{diag}\{M_1,M_2,...,M_r\},{\∀M}_i>0,i\∈[1,r]$ 表示两个收敛因子矩阵即自适用律，能够影响故障估计误差的收敛速度；P是一个对称、正定的常矩阵。

本发明提供的民机飞控系统抗干扰时变故障诊断方法与现有技术相比具有以下优点：

①本方法不是基于理想解析模型的，换言之，在设计阶段，外部干扰已经被事先考虑进来，因此提高了后续故障诊断的准确率。②可以同时获得系统状态变量和故障的具体形态，而不是仅仅给出超限警告。③本方法对故障类型不敏感，即对故障类型无确定性要求。④外部干扰和系统内部组件故障均可以是快速时变的(非常值或慢时变类型)。⑤本方法是基于非线性系统模型的，而不是传统的基于线性模型的故障诊断方法。⑥算 法规模小、运算简单，易于工程实现。本技术可进一步推广应用于一般非线性控制系统的故障诊断。

附图说明

图1为本发明提供的民机飞控系统抗干扰时变故障诊断方法流程图。

图2(a)为实验1中对状态变量x₁的估计效果图。

图2(b)为实验1中对状态变量x₁的估计效果图(局部放大图)。

图3(a)为实验1中对状态变量x₂的估计效果图。

图3(b)为实验1中对状态变量x₂的估计效果图(局部放大图)。

图4为实验1中1号故障的估计效果图。

图5为实验1中2号故障的估计效果图。

图6(a)为实验2中收敛因子矩阵取Γ₁,Μ₁时对状态变量x₁的估计效果图。

图6(b)为实验2中收敛因子矩阵取Γ₁,Μ₁时对状态变量x₁的估计效果图(放大图)。

图7(a)为实验2中收敛因子矩阵取Γ₁,Μ₁时对状态变量x₂的估计效果图。

图7(b)为实验2中收敛因子矩阵取Γ₁,Μ₁时对状态变量x₂的估计效果图(放大图)。

图8(a)为实验2中收敛因子矩阵取Γ₁,Μ₁时对1号故障的估计效果图。

图8(b)为实验2中收敛因子矩阵取Γ₁,Μ₁时对1号故障的估计效果图(放大图)。

图9(a)为实验2中收敛因子矩阵取Γ₁,Μ₁时对2号故障的估计效果图。

图9(b)为实验2中收敛因子矩阵取Γ₁,Μ₁时对2号故障的估计效果图(放大图)。

图10(a)为实验2中收敛因子矩阵取Γ₂,Μ₂时对状态变量x₁的估计效果图。

图10(b)为实验2中收敛因子矩阵取Γ₂,Μ₂时对状态变量x₁的估计效果图(放大图)。

图11(a)为实验2中收敛因子矩阵取Γ₂,Μ₂时对状态变量x₂的估计效果图。

图11(b)为实验2中收敛因子矩阵取Γ₂,Μ₂时对状态变量x₂的估计效果图(放大图)。

图12(a)为实验2中收敛因子矩阵取Γ₂,Μ₂时对1号故障的估计效果图。

图12(b)为实验2中收敛因子矩阵取Γ₂,Μ₂时对1号故障的估计效果图(放大图)。

图13(a)为实验2中收敛因子矩阵取Γ₂,Μ₂时对2号故障的估计效果图。

图13(b)为实验2中收敛因子矩阵取Γ₂,Μ₂时对2号故障的估计效果图(放大图)。

图14为实验3中收敛因子矩阵取Γ₁,Μ₁时对状态变量x₁的估计效果图。

图15为实验3中收敛因子矩阵取Γ₁,Μ₁时对状态变量x₂的估计效果图。

图16(a)为实验3中收敛因子矩阵取Γ₁,Μ₁时对1号故障的估计效果图。

图16(b)为实验3中收敛因子矩阵取Γ₁,Μ₁时对1号故障的估计效果图(放大图)。

图17(a)为实验3中收敛因子矩阵取Γ₁,Μ₁时对2号故障的估计效果图。

图17(b)为实验3中收敛因子矩阵取Γ₁,Μ₁时对2号故障的估计效果图(放大图)。

具体实施方式

下面结合附图和具体实施例对本发明提供的民机飞控系统抗干扰时变故障诊断方法进行详细说明。

如图1所示，本发明提供的民机飞控系统抗干扰时变故障诊断方法包括按顺序进行的下列步骤：

(1)建立同时包含外部干扰和内部故障向量的民机飞控系统非线性数学模型；

所述的民机飞控系统非线性数学模型是：

$\left\{\begin{array}{c}\stackrel{.}{x}\left(t\right)=\mathrm{\φ}(x,t)+B[u\left(t\right)+\mathrm{\γ}\left(t\right)]+\mathrm{Ed}\left(t\right)\\ y\left(t\right)=\mathrm{Cx}\left(t\right)\end{array}\right.---\left(1\right)$

其中，表示状态向量；表示控制输入；表示输出向量；表示系统受到的外部时变干扰；φ(x,t)是一个关于x(t)的非线性函数；表示未知的待诊断时变故障；B，C和E均为适当维的常值矩阵。

(2)设定故障估计算法的约束条件；

所述的故障估计算法的约束条件是：

①对于状态向量x(t)，非线性函数φ(x,t)满足Lipschitz条件，即存在一个常数δ₁>0使得如下不等式成立：

||φ(x₁)-φ(x₂)||≤δ₁||x₁-x₂||   (2)

注：为了简化书写，将x(t)缩写成x。此写法同样适用于后文。

②外部时变干扰d的能量是有界的，即：

||Ed||≤η₀   (3)

其中，η₀≥0表示系统能够承受的最大干扰量。

③未知的待诊断的时变故障向量γ的能量和变化频率均有界，即：

||γ||≤η₁

                                               (4) 

$\left|\right|\stackrel{.}{\mathrm{\γ}}\left|\right|\≤{\mathrm{\η}}_2$

其中，η₁≥0和η₂≥0分别表示发生故障的最大幅值和最快频率。

④状态向量x的能量和变化频率均有界，即：

||x||≤η₃

                                               (5) 

$\left|\right|\stackrel{.}{x}\left|\right|\≤{\mathrm{\η}}_4$

其中，η₃≥0和η₄≥0分别表示状态变量的最大幅值和最大变化速率。

(3)设计状态观测器；

所述的状态观测器的设计方法是：

利用步骤(1)得到的民机飞控系统非线性数学模型和步骤(2)得到的四个约束条件，设计如下式(6)所示的状态观测器，用于完成外部干扰和内部时变故障同时存在时系统状态变量_x的估计。

$\left\{\begin{array}{c}\stackrel{\stackrel{.}{^}}{x}=\mathrm{\φ}\left(\hat{x}\right)+B[u+\widehat{\mathrm{\γ}}]+L(y-\hat{y})\\ \hat{y}=C\hat{x}\end{array}\right.---\left(6\right)$

其中，表示对系统状态变量x的估计值；表示对时变故障γ的估计值； 表示对系统输出y的估计值；L是一个可调节增益矩阵。

(4)设计自适应故障估计算法；

所述的自适应故障估计算法的设计方法是：

结合步骤(1)得到的民机飞控系统非线性数学模型和步骤(3)得到的状态观测器，设计如下式(7)所示的自适应故障估计算法，用于完成对系统内部时变故障γ的精确估计。

$\stackrel{\stackrel{.}{^}}{\mathrm{\γ}}={\mathrm{\ΓB}}^{T}P\tilde{x}+M\stackrel{\stackrel{.}{~}}{x}---\left(7\right)$

其中，表示状态观测器的观测误差；表示误差估计的误差；  $\mathrm{\Γ}=\mathrm{diag}\{{\mathrm{\Γ}}_1,{\mathrm{\Γ}}_2,...,{\mathrm{\Γ}}_r\},\∀{\mathrm{\Γ}}_i>0,i\∈[1,r]$ 和 $M=\mathrm{diag}\{M_1,M_2,...,M_r\},{\∀M}_i>0,i\∈[1,r]$ 表示两个收敛因子矩阵即自适用律，能够影响故障估计误差的收敛速度；P是一个对称、正定的常矩阵。

为了验证本发明提供的民机飞控系统抗干扰时变故障诊断方法的稳定性，本发明人 对其进行了证明，过程如下：

定理1如果式(7)中的P矩阵和对称正定矩阵Q满足以下条件：

(-LC)^TP+P(-LC)+δPP+δI＝-Q   (8)

其中δ>0是满足Lipschitz条件的正数，I是适当维度的单位阵。那么，基于状态观测器和自适应故障估计算法的民机飞控系统抗干扰时变故障诊断方法能够保证状态估计误差和误差估计误差收敛于有限的正数，即：

$\underset{t\→\∞}{\lim }\left|\right|\tilde{x}\left(t\right)\left|\right|\≤{\mathrm{\ξ}}_1---\left(9\right)$

$\underset{t\→\∞}{\lim }\left|\right|\widetilde{\mathrm{\γ}}\left(t\right)\left|\right|\≤{\mathrm{\ξ}}_2---\left(10\right)$

其中ξ₁>0，ξ₂>0是两个有限的正数。

证明：定义一个Lyapunov函数V(t)如下：

$V\left(t\right)={\tilde{x}}^{T}P\tilde{x}+{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\widetilde{\mathrm{\γ}}---\left(11\right)$

注：V(t)能够同时衡量状态估计误差和误差估计误差

将式(11)对时间求导，得到下面的表达式：

$\stackrel{.}{V}={\stackrel{\stackrel{.}{~}}{x}}^{T}P\tilde{x}+{\tilde{x}}^{T}P\stackrel{\stackrel{.}{~}}{x}+{\stackrel{\stackrel{.}{~}}{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\widetilde{\mathrm{\γ}}+{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\stackrel{\stackrel{.}{~}}{\mathrm{\γ}}---\left(12\right)$

在中共有4项，为了便于后续阅读与分析，将其分别展开成如下形式：

${\stackrel{\stackrel{.}{~}}{x}}^{T}P\tilde{x}=[{\tilde{x}}^T{(-\mathrm{LC})}^T+{(\mathrm{\φ}-\widehat{\mathrm{\φ}})}^T+{\widetilde{\mathrm{\γ}}}^T{B}^T+d^T{E}^T]P\tilde{x}---\left(13\right)$

${\tilde{x}}^{T}P\stackrel{\stackrel{.}{~}}{x}={\tilde{x}}^{T}P[(-\mathrm{LC})\tilde{x}+(\mathrm{\φ}-\widehat{\mathrm{\φ}})+B\widetilde{\mathrm{\γ}}+\mathrm{Ed}]---\left(14\right)$

$\begin{array}{c}{\stackrel{\stackrel{.}{~}}{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\widetilde{\mathrm{\γ}}={(\stackrel{.}{\mathrm{\γ}}-\stackrel{\stackrel{.}{^}}{\mathrm{\γ}})}^T{\mathrm{\Γ}}^{-1}\widetilde{\mathrm{\γ}}\\ ={\stackrel{.}{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\widetilde{\mathrm{\γ}}+{\left({-\mathrm{\ΓB}}^{T}P\tilde{x}\right)}^T{\mathrm{\Γ}}^{-1}\widetilde{\mathrm{\γ}}+{(-M\stackrel{\stackrel{.}{~}}{x})}^T{\mathrm{\Γ}}^{-1}\widetilde{\mathrm{\γ}}\\ ={\stackrel{.}{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\widetilde{\mathrm{\γ}}-{\tilde{x}}^T\mathrm{PB}\widetilde{\mathrm{\γ}}-{\stackrel{\stackrel{.}{~}}{x}}^{T}M{\mathrm{\Γ}}^{-1}\widetilde{\mathrm{\γ}}\end{array}---\left(15\right)$

$\begin{array}{c}{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\stackrel{\stackrel{.}{~}}{\mathrm{\γ}}={\left({\stackrel{\stackrel{.}{~}}{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\widetilde{\mathrm{\γ}}\right)}^T\\ ={\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\stackrel{.}{\mathrm{\γ}}-{\widetilde{\mathrm{\γ}}}^T{B}^{T}P\tilde{x}-{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}M\stackrel{\stackrel{.}{~}}{x}\end{array}---\left(16\right)$

将式(13)(14)(15)(16)代入式(12)后，变为下式：

$\begin{array}{c}\stackrel{.}{V}={\tilde{x}}^T[{(-\mathrm{LC})}^{T}P+P(-\mathrm{LC})]\tilde{x}+[{(\mathrm{\φ}-\widehat{\mathrm{\φ}})}^{T}P\tilde{x}+{\tilde{x}}^{T}P(\mathrm{\φ}-\widehat{\mathrm{\φ}})]+\\ [d^T{E}^{T}P\tilde{x}+{\tilde{x}}^T\mathrm{PEd}]+[2\widetilde{\mathrm{\γ}}{\mathrm{\Γ}}^{-1}\stackrel{.}{\mathrm{\γ}}-2{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}M\stackrel{\stackrel{.}{~}}{x}]\end{array}---\left(17\right)$

在应用Lipschitz条件之后，可以得到如下不等式：

$\begin{array}{c}\stackrel{.}{V}\≤{\tilde{x}}^T[{(-\mathrm{LC})}^{T}P+P(-\mathrm{LC})]\tilde{x}+2\left|\right|{\tilde{x}}^{T}P\left|\right|\mathrm{\δ}\left|\right|\tilde{x}\left|\right|+\\ 2\left|\right|\mathrm{Ed}\left|\right|{\mathrm{\λ}}_{\max }\left(P\right)\left|\right|\tilde{x}\left|\right|+[2{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\stackrel{.}{\mathrm{\γ}}-2{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}M\stackrel{\stackrel{.}{~}}{x}]\end{array}---\left(18\right)$

其中λ_max(P)是P矩阵的最大特征值。

接下来，将式(18)中的最后一项整理为如下形式：

$2{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\stackrel{.}{\mathrm{\γ}}-2{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}M\stackrel{\stackrel{.}{~}}{x}=2{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}(\stackrel{.}{\mathrm{\γ}}-M\stackrel{\stackrel{.}{~}}{x})---\left(19\right)$

这里，定义为如下形式： 

$\stackrel{.}{Z}=\stackrel{.}{\mathrm{\γ}}-M\stackrel{\stackrel{.}{~}}{x}---\left(20\right)$

将式(20)代入式(19)后，得到： 

$\begin{array}{c}2{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\stackrel{.}{\mathrm{\γ}}-2{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}M\stackrel{\stackrel{.}{~}}{x}=2{\widetilde{\mathrm{\γ}}}^T{\mathrm{\Γ}}^{-1}\stackrel{.}{Z}\\ \≤{\widetilde{\mathrm{\γ}}}^T\widetilde{\mathrm{\γ}}+{\stackrel{.}{Z}}^T{\mathrm{\Γ}}^{-1}{\mathrm{\Γ}}^{-1}\stackrel{.}{Z}\\ \≤{\left|\right|\widetilde{\mathrm{\γ}}\left|\right|}^2+{\left|\right|\stackrel{.}{Z}\left|\right|}^2{\mathrm{\λ}}_{\max }\left({\mathrm{\Γ}}^{-1}{\mathrm{\Γ}}^{-1}\right)\end{array}---\left(21\right)$

很明显，将式(21)代入式(18)之后，式(18)变为下式： 

$\begin{array}{c}\stackrel{.}{V}\≤{\tilde{x}}^T[{(-\mathrm{LC})}^{T}P+P(-\mathrm{LC})]\tilde{x}+\mathrm{\δ}[{\left|\right|{\tilde{x}}^{T}P\left|\right|}^2+{\left|\right|\tilde{x}\left|\right|}^2]+\\ 2{\mathrm{\η}}_0{\mathrm{\λ}}_{\max }\left(P\right)\left|\right|\tilde{x}\left|\right|+{\left|\right|\widetilde{\mathrm{\γ}}\left|\right|}^2+{\left|\right|\stackrel{.}{Z}\left|\right|}^2{\mathrm{\λ}}_{\max }\left({\mathrm{\Γ}}^{-1}{\mathrm{\Γ}}^{-1}\right)\\ ={\tilde{x}}^T[{(-\mathrm{LC})}^{T}P+P(-\mathrm{LC})]\tilde{x}+\mathrm{\δ}[{\tilde{x}}^T\mathrm{PP}\tilde{x}+{\tilde{x}}^T\tilde{x}]+\\ 2{\mathrm{\η}}_0{\mathrm{\λ}}_{\max }\left(P\right)\left|\right|\tilde{x}\left|\right|+{\left|\right|\widetilde{\mathrm{\γ}}\left|\right|}^2+{\left|\right|\stackrel{.}{Z}\left|\right|}^2{\mathrm{\λ}}_{\max }\left({\mathrm{\Γ}}^{-1}{\mathrm{\Γ}}^{-1}\right)\\ ={\tilde{x}}^T[{(-\mathrm{LC})}^{T}P+P(-\mathrm{LC})+\mathrm{\δPP}+\mathrm{\δI}]\tilde{x}+\\ 2{\mathrm{\η}}_0{\mathrm{\λ}}_{\max }\left(P\right)\left|\right|\tilde{x}\left|\right|+{\left|\right|\widetilde{\mathrm{\γ}}\left|\right|}^2+{\left|\right|\stackrel{.}{Z}\left|\right|}^2{\mathrm{\λ}}_{\max }\left({\mathrm{\Γ}}^{-1}{\mathrm{\Γ}}^{-1}\right)\end{array}---\left(22\right)$

其中η₀参见式(3)的定义。

因为||γ||≤η₁，则故障估计值的分布区间是[-η₁,+η₁]，可以得到下述不等式：

$\left|\right|\widetilde{\mathrm{\γ}}\left|\right|\≤{2\mathrm{\η}}_1---\left(23\right)$

其中η₁参见式(4)的定义。

另外，由于和均有界，也一定有界。因此，定义：

$\left|\right|\stackrel{.}{Z}\left|\right|\≤{\mathrm{\η}}_5---\left(24\right)$

其中，η₅可以取η₅＝η₂+2λ_max(Μ)η₄

应用式(23)和式(24)之后，式(22)可以进一步变换为下式：

$\stackrel{.}{V}\≤{\tilde{x}}^T[{(-\mathrm{LC})}^{T}P+P(-\mathrm{LC})+\mathrm{\δPP}+\mathrm{\δI}]\tilde{x}+2{\mathrm{\η}}_0{\mathrm{\λ}}_{\max }\left(P\right)\left|\right|\tilde{x}\left|\right|+4{\mathrm{\η}}_1^2+{\mathrm{\η}}_5^2{\mathrm{\λ}}_{\max }\left({\mathrm{\Γ}}^{-1}{\mathrm{\Γ}}^{-1}\right)---\left(25\right)$

为了便于阅读，定义常数C如下：

$C=4{\mathrm{\η}}_1^2+{\mathrm{\η}}_5^2{\mathrm{\λ}}_{\max }\left({\mathrm{\Γ}}^{-1}{\mathrm{\Γ}}^{-1}\right)---\left(\text{26}\right)$

将式(8)和式(26)代入式(25)之后，不等式被化为如下形式：

$\begin{array}{c}\stackrel{.}{V}\≤-{\tilde{x}}^{T}Q\tilde{x}+2{\mathrm{\η}}_0{\mathrm{\λ}}_{\max }\left(P\right)\left|\right|\tilde{x}\left|\right|+C\\ \≤-{\mathrm{\λ}}_{\min }\left(Q\right){\left|\right|\tilde{x}\left|\right|}^2+2{\mathrm{\η}}_0{\mathrm{\λ}}_{\max }\left(P\right)\left|\right|\tilde{x}\left|\right|+C\end{array}---\left(27\right)$

其中λ_min(Q)代表满足定理1的正定矩阵Q的最小特征值。

接下来，为了完成稳定性分析，分两种情况进行讨论。

情况1： $\left|\right|\tilde{x}\left|\right|>\mathrm{\σ}$ 其中 $\mathrm{\σ}=\frac{{2\mathrm{\η}}_0{\mathrm{\λ}}_{\max }\left(P\right)+\sqrt{4{{\mathrm{\η}}_0}^2{\mathrm{\λ}}_{\max }^2\left(P\right)+4C{\mathrm{\λ}}_{\min }\left(Q\right)}}{2{\mathrm{\λ}}_{\min }\left(\text{Q}\right)}$

在这种情况下，可以得出从而确定公式(6)所描述的就是非线性故障系统(参见公式(1)的描述)的稳定状态观测器。

情况2： $\left|\right|\tilde{x}\left|\right|>\mathrm{\σ}$

在这种情况下，可以得出且状态观测器是发散的，这会造成增大，使得 从而又满足情况1，情况1又会使得减小，以此类推，也就是说状态误差会稳定在一个确定的误差范围之内。

综合考虑情况1和情况2，状态观测器是稳定的，并且所设计的状态观测器能够以一定精确度跟踪状态变量x(t)和时变故障γ(t)。

注：式(8)中的L是状态观测器的收敛矩阵，L只要能够保证(-LC)是稳定矩阵就可以。

例如，下面的L矩阵是合理的：

L＝PC^T   (28)

另外，本发明人还对本发明提供的民机飞控系统抗干扰时变故障诊断方法进行了实验，以对其效果进行的评估，过程如下：

为了充分验证本发明提供的民机飞控系统抗干扰时变故障诊断方法的有效性，选取具有如下参数形式的非线性模型：

$\mathrm{\φ}\left(x\right)=\left[\begin{array}{c}{x}_2\\ -0.5\sin \left(x_1\right)\end{array}\right]---\left(29\right)$

$B=\left[\begin{array}{cc}1.5& 0\\ 0& 1\end{array}\right],E=C=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$

显然，φ(x)满足Lipschitz条件。根据式(8)建立具有如下参数形式的状态观测器：

δ＝1.0

$P=L=\left[\begin{array}{cc}4.5& -3.5\\ -3.5& 5.5\end{array}\right],Q=\left[\begin{array}{cc}31.5& -35\\ -35& 41.5\end{array}\right]---\left(\text{30}\right)$

实验1不存在内部故障和外部干扰的实验

假设在故障诊断过程中没有故障和外部干扰，即

选取具有如下形式的控制输入向量u：

$u=\left[\begin{array}{c}0.1\\ 0.5\end{array}\right]---\left(31\right)$

注：通过实验我们发现：输入向量u的形式变化对于状态和故障估计没有显著影响，所以，我们这里选取了一个简单的常数向量。

初始状态向量x和估计向量的选择如下所示：

$x=\left[\begin{array}{c}10.0\\ 20.0\end{array}\right],\hat{x}=\left[\begin{array}{c}0.0\\ 0.0\end{array}\right]---\left(32\right)$

自适应故障估计算法中的收敛因子矩阵的选择如下所示：

${\mathrm{\Γ}}_1=\left[\begin{array}{cc}0.5& 0\\ 0& 0.5\end{array}\right],M_1=\left[\begin{array}{cc}5.0& 0\\ 0& 18.0\end{array}\right]---\left(\text{33}\right)$

仿真过程持续10秒钟。状态估计效果参见图2(a)、2(b)、3(a)、3(b)。故障估计效果参见图4、5。

实验2：存在外部干扰和内部故障的实验

假设在故障诊断过程中同时存在式(34)形式的外部时变干扰和式(35)形式的内部时变故障。控制输入向量仍然选用式(31)所给出的形式，初始状态向量x和估计向量^x选用式(32)所给出的形式。

$d=\left[\begin{array}{c}0.1\cos (62.8\·t)\\ 0.5\sin (62.8\·t)\end{array}\right]---\left(34\right)$

$\mathrm{\γ}=\left[\begin{array}{c}0.08\cos (6.28\·t)\\ 0.4\sin (6.28\·t)\end{array}\right]---\left(35\right)$

为了充分分析和评价状态估计及故障估计的快速性，将自适应故障估计算法中的收敛因子矩阵分别选为式(33)和式(36)两种形式。

${\mathrm{\Γ}}_2=\left[\begin{array}{cc}0.1& 0\\ 0& 0.1\end{array}\right],M_2=\left[\begin{array}{cc}3.0& 0\\ 0& 9.0\end{array}\right]---\left(\text{36}\right)$

基于Γ₁和Μ₁的估计效果参见图6(a)、6(b)、7(a)、7(b)、8(a)、8(b)、9(a)、9(b)。

基于Γ₂和Μ₂的估计效果参见图10(a)、10(b)、11(a)、11(b)、12(a)、12(b)、13(a)、13(b)。

实验3：存在系统内部故障，不存在外部干扰的实验

假设在故障诊断过程中不存在外部干扰，但存在式(35)中描述的系统内部故障。控制输入向量与式(31)定义的输入向量相同；初始状态向量x和估计向量与式(32)定义的相同。

仿真过程同样持续10秒钟。自适应故障估计算法中的收敛因子按照式(33)选取。状态估计效果参见图14、15。故障估计效果参见图16(a)、16(b)、17(a)、17(b)。

综合实验1到实验3的实验结果，可以得出如下结论：

1)状态估计的效果很理想。不管外部干扰或者内部故障存在与否，状态估计的误差稳态值都非常小，即另外，状态估计的快速性也非常好，渐近到x的过程耗时在0.4秒以内(参见图2，图3，图6，图7，图14，图15)。这表明状态估计的结果是令人满意的；

2)外部干扰对于故障估计效果有明显的影响。通过对比实验2和实验3，可以发现：当存在外部时变干扰时，在线故障估计值有明显的小幅振荡(参见图8(b)，图9(b)，图12(b)和图13(b))。在这种情况下，故障估计值会围绕在真实值附近振荡。另外，振荡均值所代 表的正是实际的故障值。当没有外部干扰时，所估计的故障值会无震荡地收敛到实际值。(参见图16(b)和图17(b))。

3)收敛因子对故障估计的快速性有显著影响。在实验2中，通过比较Γ₁,Μ₁和Γ₂,Μ₂，可以发现：如果收敛因子取Γ₁,Μ₁，故障估计值会在大约6秒内进入到稳定振荡阶段(参见图8(b))。然而，当收敛因子选为Γ₂,Μ₂时，这个过程会被延长至大约15秒(参见图12(b))。因此，适当地增大收敛因子可以提高故障诊断的快速性。然而，同样通过实验还发现：太大的收敛因子会导致状态观测器发散或非预期的超调量，因此，收敛因子必须谨慎选择。

Claims (5)

1.一种民机飞控系统抗干扰时变故障诊断方法，其特征在于：其包括按顺序进行的下列步骤：

(1)建立同时包含外部干扰和内部故障向量的民机飞控系统非线性数学模型；

(2)设定故障估计算法的约束条件；

(3)设计状态观测器；

(4)设计自适应故障估计算法。

2.根据权利要求1所述的民机飞控系统抗干扰时变故障诊断方法，其特征在于：在步骤(1)中，所述的民机飞控系统非线性数学模型是：

$\left\{\begin{array}{c}\stackrel{\·}{x}\left(t\right)=\mathrm{\φ}(x,t)+B[u\left(t\right)+\mathrm{\γ}\left(t\right)]+\mathrm{Ed}\left(t\right)\\ y\left(t\right)=\mathrm{Cx}\left(t\right)\end{array}\right.---\left(1\right)$

其中，表示状态向量；表示控制输入；表示输出向量；表示系统受到的外部时变干扰；φ(x,t)是一个关于x(t)的非线性函数；表示未知的待诊断时变故障；B，C和E均为适当维的常值矩阵。

3.根据权利要求1所述的民机飞控系统抗干扰时变故障诊断方法，其特征在于：在步骤(2)中，所述的故障估计算法的约束条件是：

①对于状态向量x(t)，非线性函数φ(x,t)满足Lipschitz条件，即存在一个常数δ₁>0使得如下不等式成立：

||φ(x₁)-φ(x₂)||≤δ₁||x₁-x₂||   (2)

x(t)缩写成x；

②外部时变干扰d的能量是有界的，即：

||Ed||≤η₀   (3)

其中，η₀≥0表示系统能够承受的最大干扰量；

③未知的待诊断的时变故障向量γ的能量和变化频率均有界，即：

$\begin{array}{c}\left|\right|\mathrm{\γ}\left|\right|\≤{\mathrm{\η}}_1\\ \left|\right|\stackrel{\·}{\mathrm{\γ}}\left|\right|\≤{\mathrm{\η}}_2\end{array}---\left(4\right)$

其中，η₁≥0和η₂≥0分别表示发生故障的最大幅值和最快频率；

④状态向量x的能量和变化频率均有界，即：

$\begin{array}{c}\left|\right|x\left|\right|\≤{\mathrm{\η}}_3\\ \left|\right|\stackrel{\·}{x}\left|\right|\≤{\mathrm{\η}}_4\end{array}---\left(5\right)$

其中，η₃≥0和η₄≥0分别表示状态变量的最大幅值和最大变化速率。

4.根据权利要求1所述的民机飞控系统抗干扰时变故障诊断方法，其特征在于：在步骤(3)中，所述的状态观测器的设计方法是：

利用步骤(1)得到的民机飞控系统非线性数学模型和步骤(2)得到的四个约束条件，设计如下式(6)所示的状态观测器，用于完成外部干扰和内部时变故障同时存在时系统状态变量_x的估计；

$\left\{\begin{array}{c}\stackrel{\stackrel{\·}{^}}{x}=\mathrm{\φ}\left(\hat{x}\right)+B[u+\widehat{\mathrm{\γ}}]+L(y-\hat{y})\\ \hat{y}=C\hat{x}\end{array}\right.---\left(6\right)$

其中，表示对系统状态变量x的估计值；表示对时变故障γ的估计值；

表示对系统输出y的估计值；L是一个可调节增益矩阵。

5.根据权利要求1所述的民机飞控系统抗干扰时变故障诊断方法，其特征在于：在步骤(4)中，所述的自适应故障估计算法的设计方法是：

结合步骤(1)得到的民机飞控系统非线性数学模型和步骤(3)得到的状态观测器，设计如下式(7)所示的自适应故障估计算法，用于完成对系统内部时变故障γ的精确估计：

$\stackrel{\stackrel{\·}{^}}{\mathrm{\γ}}=\mathrm{\Γ}{B}^{T}P\tilde{x}+M\stackrel{\stackrel{\·}{~}}{x}---\left(7\right)$

其中，表示状态观测器的观测误差；表示误差估计的误差； $\mathrm{\Γ}=\mathrm{diag}\{{\mathrm{\Γ}}_1,{\mathrm{\Γ}}_2,...,{\mathrm{\Γ}}_r\}\∀{\mathrm{\Γ}}_i>0,i\∈[1,r]$ 和 $M=\mathrm{diag}\{M_1,M_2,...,M_r\}\∀M_i>0,i\∈[1,r]$ 表示两个收敛因子矩阵即自适用律，能够影响故障估计误差的收敛速度；P是一个对称、正定的常矩阵。