CN101807047B - 基于模糊奇偶方程及ar模型的非线性系统故障预测方法 - Google Patents

Abstract

一种基于模糊奇偶方程及AR模型的非线性系统故障预测方法，包括以下步骤：采用模糊奇偶方程方法估计非线性系统执行器或传感器偏差；采用AR模型对所述执行器、传感器产生的偏差序列进行建模，给出偏差预测值；由偏差预测值结合其统计规律计算执行器或传感器故障发生概率，并用预测置信因子对预测准确性进行评价。从而提供一种以故障发生概率的形式给出了预测结果，并预测置信因子来反映了预测结果的准确程度的基于模糊奇偶方程及AR模型的非线性系统故障预测方法。

Description

基于模糊奇偶方程及AR模型的非线性系统故障预测方法

技术领域

本发明涉及一种故障诊断与预测领域，特别是一种基于模糊奇偶方程及AR模型的非线性系统故障预测方法。

背景技术

奇偶方程是线性系统故障诊断中常用的一种基于模型的方法，主要是利用实测数据去检测系统数学方程的一致性。奇偶方程产生的残差在理论上仅在故障时非零，由此可进行故障诊断。将T-S模糊模型、全解耦奇偶方程和参数估计相结合，同时对非线性系统多个传感器故障进行检测、隔离和识别，通过对飞机控制系统传感器的故障诊断验证了方法的有效性。时间序列方法首先由Box和Jenkens在1970年系统提出，发展到现在已较为成熟，在社会学、自然科学以及工程等众多学科领域都取得了良好的效果，能解决社会科学中诸如供应链管理(SCM)中的需求预报、不同季节的大气污染指标(API)预报等问题；在对设备的故障监测、诊断和预报中的运用也极其广泛。Fassois S.D.和Sakellariou J.S.运用时间序列分析方法对飞机仪表板的故障情况进行识别；吴庚申等人采用Bently实验台所采集的碰摩、松动、不对中和不平衡四种典型汽轮机转子振动故障水平方向与垂直方向的数据，建立了汽轮机转子振动故障序列自回归滑移平均(ARMA)模型。

而现有技术基于模糊奇偶方程的诊断方法，对于故障参数估计都针对于常值故障及刻度系数故障，对于故障幅值随时间变化的情况还较少应用，且使用时间序列分析方法进行预测在故障预测方面应用较少，对于预测结果的准确性没有给出反映指标。

发明内容

针对上述现有技术的缺陷，本发明的目的是提供一种以故障发生概率的形式给出了预测结果，并预测置信因子来反映了预测结果的准确程度的基于模糊奇偶方程及AR模型的非线性系统故障预测方法。

为达到上述目的，本发明采用如下技术方案：

一种基于模糊奇偶方程及AR模型的非线性系统故障预测方法，包括以下步骤：

采用模糊奇偶方程方法估计非线性系统执行器或传感器偏差；

采用AR模型对所述执行器、传感器产生的偏差序列进行建模，给出偏差预测值；

由偏差预测值结合其统计规律计算执行器或传感器故障发生概率，并用预测置信因子对预测准确性进行评价。

本发明的基于模糊奇偶方程及AR模型的非线性系统故障预测方法，其中所述采用模糊奇偶方程方法估计非线性系统执行器、传感器偏差的步骤还包括：

对非线性系统在各工作点处线性化，得到局部线性模型；

离线计算出各局部线性模型的奇偶方程；

利用T-S模型产生系统的残差，根据残差进行故障诊断。

本发明的基于模糊奇偶方程及AR模型的非线性系统故障预测方法，其中所述利用T-S模型产生系统的残差，根据残差进行故障诊断步骤还包括：执行器残差的计算和故障征兆评估，以及传感器残差的计算和故障征兆评估。

本发明的基于模糊奇偶方程及AR模型的非线性系统故障预测方法，其中所述采用AR模型对所述执行器、传感器产生的偏差序列进行建模的步骤中，采用滚动数据窗方法给残差建模，并采用FPE准则进行模型定阶。

本发明的基于模糊奇偶方程及AR模型的非线性系统故障预测方法，其中所述采用AR模型对所述执行器、传感器产生的偏差序列进行建模的步骤中还包括对模型进行平稳性检验，对模型残差进行白噪声检验的步骤。

本发明的基于模糊奇偶方程及AR模型的非线性系统故障预测方法，其中所述置信因子用于反映由于预测步长增加等因素而导致的预测准确度的降低程度。

本发明的基于模糊奇偶方程及AR模型的非线性系统故障预测方法，其中所述故障预测包括预测正确率和故障正确预测率两个主要评价参数。

本发明的基于模糊奇偶方程及AR模型的非线性系统故障预测方法，其中所述预测置信因子用c_k表示，且

其中ξ_e(t+k)＝3σ_et(1)，σ_et(k)为k步的预测误差。

本发明的基于模糊奇偶方程及AR模型的非线性系统故障预测方法，其中所述ξ_e(t+k)设为ξ_e(t+k)＝|x(t+k)-T_D|，其中x(t+k)为t+k时刻残差值，T_D为阈值。

本发明基于现有故障预测方法的不足，针对工程实际中普遍存在的非线性系统，给出了基于模糊奇偶方程和模糊基函数网络的故障预测方法，并针对预测的不准确性，以故障发生概率的形式给出了预测结果，能够给人以更直观的印象，并以预测置信因子来反映了预测结果的准确程度。

附图说明

图1是本发明基于模糊奇偶方程及AR模型的非线性系统故障预测方法的方法流程图。

具体实施方式

下面结合附图对本发明基于模糊奇偶方程及AR模型的非线性系统故障预测方法的实施方式进行详细说明。

参见图1，步骤S1，采用模糊奇偶方程方法估计非线性系统执行器或传感器偏差。具体实现方法如下：

考虑非线性系统

$\left\{\begin{array}{c}\stackrel{\·}{x}\left(t\right)=f[x\left(t\right),u\left(t\right),w\left(t\right)]\\ y\left(t\right)=h[x\left(t\right),u\left(t\right),w\left(t\right)]\end{array}\right.---\left(1\right)$

x(t)∈Rⁿ为系统状态，u(t)∈R^p为执行器输入，y(t)∈R^q为传感器输出，w(t)∈R^r为扰动输入，f[x(t)，u(t)，w(t)和h[x(t)，u(t)，w(t)]为光滑的非线性函数。则可利用T-S模糊模型将此非线性系统在一系列工作点附近线性化得到局部线性模型，每一个局部线性模型都可以描述对应工作点处局部系统性能，全局系统性能由所有局部线性模型的输出融合而成。

设非线性系统(1)在工作点l(l＝1，2，…，m)处线性化，可得：

$\left\{\begin{array}{c}x(k+1)=A_{l}x\left(k\right)+B_{l}u\left(k\right)+F_{l}w\left(k\right)\\ y\left(k\right)=C_{l}x\left(k\right)+D_{l}u\left(k\right)+G_{l}w\left(k\right)\end{array}\right.---\left(2\right)$

式中l对应系统的工作点l。x(k)，u(k)，y(k)，w(k)系统的状态、执行器输入、传感器输出和扰动输入。A_l、B_l、C_l、D_l、F_l和G_l是合适维数的矩阵。

对于局部线性模型(2)建立奇偶方程。设系统(1)在工作点l处的奇偶方程残差为r_l(k)(l＝1，2，…，m)，则可用T-S模型的“如果…则…”(IF-THEN)模糊规则描述如下：

规则l(l＝1，2，…，m)：

如果：ψ₁(k)为S_L1且ψ₂(k)为S_l2且…，则系统(1)的奇偶方程残差为r_l(k)。

其中m为局部线性模型的个数(规则个数)，每一条规则对应一个工作点。ψ＝[ψ₁ψ₂…ψ_ρ]^T为T-S模型的前件变量，S_lj(j＝1，2，…，ρ)为模糊集合。前件变量ψ_j(k)对模糊集S_lj的隶属函数(隶属度)为

系统残差即为上述T-S模型的输出：

$r\left(k\right)=\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l^{*}\left(\mathrm{\ψ}\left(k\right)\right)\·r_l\left(k\right)---\left(3\right)$

式中β_l ^*(ψ(k))需满足：

$\left\{\begin{array}{c}\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l^{*}\left(\mathrm{\ψ}\left(k\right)\right)=1\\ 0\≤{\mathrm{\β}}_l^{*}\left(\mathrm{\ψ}\left(k\right)\right)\≤1\end{array}\right.---\left(4\right)$

可令

${\mathrm{\β}}_l^{*}\left(\mathrm{\ψ}\left(k\right)\right)=\frac{1}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)---\left(5\right)$

其中β_l(ψ(k))为模糊规则l的执行度，可用如下公式计算：

${\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)=\underset{j=1}{\overset{\mathrm{\ρ}}{\mathrm{\Π}}}{\mathrm{\μ}}_{s_{\mathrm{lj}}}\left({\mathrm{\ψ}}_j\left(k\right)\right)---\left(6\right)$

则式(3)可写为

$r\left(k\right)=\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·r_l\left(k\right)}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}---\left(7\right)$

式(7)称为非线性系统(2)的模糊奇偶方程。

对于非线性系统(2)，可在各工作点处线性化得到局部线性模型，离线计算出各局部线性模型的奇偶方程，再用模糊奇偶方程(7)产生系统的残差，根据残差进行故障诊断。这样可减少在线计算量，提高计算速度。

执行器残差计算方法：

假设非线性系统中传感器工作正常，仅考虑执行器故障。

对于数据窗内s+1个最新测量数据y(k-s)至y(k)，由式(2)可得具有时间冗余的测量方程为

Y(k)＝H_l0x(k-s)+H_lcU(k)+H_lwW(k)(8)

其中，U(k)＝[u^T(k-s)…u^T(k)]^T，W(k)＝[w^T(k-s)…w^T(k)]^T，Y(k)＝[y^T(k-s)…y^T(k)]^T。下标l表示工作点。矩阵H_l0、H_lc、H_lw是(2)式所示的工作点l处的时间冗余测量矩阵。

$H_{l0}=\left[\begin{array}{c}{C}_l\\ C_l{A}_l\\ .\\ .\\ .\\ C_l{A}_l^s\end{array}\right]---\left(9\right)$

在k时刻第l个工作点对第i个执行器敏感的全解耦奇偶方程为：

$r_i^l\left(k\right)=v_{\mathrm{li}}^T[Y\left(k\right)-H_{\mathrm{lc}}{U}_c\left(k\right)]---\left(12\right)$

其中，下标i对应第i个执行器，r_i ^l(k)为对执行器i敏感的残差。U_c(k)＝[u_c(k-s)…u_c(k)]^T，u_c(k-s)，…，u_c(k)为正常的执行器输入。v_li为对第i个执行器敏感的全解耦奇偶向量。由1.2节可知，v_li应满足：

$v_{\mathrm{li}}^T\left[\begin{array}{ccc}{H}_{l0}& H_{\mathrm{lw}}& H_{\mathrm{lci}}\end{array}\right]=0---\left(13\right)$

式中H_lci即将H_lc中的D_l、B_l用D_l ^*、B_l ^*代替，D_l ^*、B_l ^*为从D_l、B_l中划去与第i个执行器对应的第i列。

假设矩阵[H_l0H_lwH_lci]中不相关的列向量数量为n_x，则式(13)有非零解的充要条件为：

$s>\frac{n_x}{q}-1---\left(14\right)$

充分条件为：

$\left\{\begin{array}{c}s>\frac{n}{q+1-(r+p)}-1\\ q>(r+p)-1\end{array}\right.---\left(15\right)$

其中n为系统状态维数、q为传感器个数、r为扰动输入的维数、p为执行器输入的维数。

式(15)说明，当传感器数(q)多于执行器数(p)与扰动输入数(r)之和时，总可找到合适的s，使式(13)的解存在。

根据式(7)可得对第i个执行器敏感的残差为

$r_i\left(k\right)=\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·\left[v_{\mathrm{li}}^T(Y\left(k\right)-H_{\mathrm{lc}}{U}_c\left(k\right))\right]}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}---\left(16\right)$

执行器故障征兆估计方法：

将式(8)代入式(7)，有：

$r_i\left(k\right)=\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·\left[v_{\mathrm{li}}^T(Y\left(k\right)-H_{\mathrm{lc}}{U}_c\left(k\right))\right]}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}---\left(17\right)$

$=\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·\left[v_{\mathrm{li}}^T(H_{l0}{x}^l(k-s)+H_{\mathrm{lc}}U\left(k\right)+H_{\mathrm{lw}}{W}_l\left(k\right)-H_{\mathrm{lc}}{U}_c\left(k\right))\right]}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}$

将式(13)代入式(17)，有：

$r_i\left(k\right)=\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·\left[v_{\mathrm{li}}^T{H}_{\mathrm{lc}}(U\left(k\right)-H_{\mathrm{lc}}{U}_c\left(k\right))\right]}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}---\left(18\right)$

$=\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{li}}^T{H}_{\mathrm{lc}}}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}U\left(k\right)-\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{li}}^T{H}_{\mathrm{lc}}}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}{U}_c\left(k\right)$

考虑执行器i的故障征兆模型，有：

u_i(k)＝η_i(k)·u_ic(k)+λ_i(k)(19)

其中u_i(k)为执行器i的实际输入，其刻度因子为η_i(k)，偏差为λ_i(k)，u_ic(k)为执行器i无故障时的输入。

设在数据窗内同一个执行器具有相同的状态，则

U(k)＝η_i(k)·U_c(k)+λ_i(k)E    (20)

式中E＝[11…1]^T，是元素全为1的s+1维向量。将(20)式代入(18)式可得：

$r_i\left(k\right)=\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{li}}^T{H}_{\mathrm{lc}}}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\mathrm{\ψ}\left(k\right)}\{{\mathrm{\η}}_i\left(k\right)\·U_c\left(k\right)+{\mathrm{\λ}}_i\left(k\right)E\}-\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{li}}^T{H}_{\mathrm{lc}}}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}{U}_c\left(k\right)$

$=({\mathrm{\η}}_i\left(k\right)-1)\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{li}}^T{H}_{\mathrm{lc}}}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}{U}_c\left(k\right)+{\mathrm{\λ}}_i\left(k\right)\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{li}}^T{H}_{\mathrm{lc}}}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}E$

$=[\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{li}}^T{H}_{\mathrm{lc}}}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}{U}_c\left(k\right),{\mathrm{\λ}}_i\left(k\right)\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{li}}^T{H}_{\mathrm{lc}}}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}E]{[({\mathrm{\η}}_i\left(k\right)-1),{\mathrm{\λ}}_i\left(k\right)]}^T$

$=\mathrm{\φ}\left(k\right)\mathrm{\θ}\left(k\right)$

考虑测量噪声，则有：

r_i(k)＝φ(k)θ(k)+n(k)(21)

其中：

θ(k)＝[(η_i(k)-1)，λ_i(k)]^T    (22)

$\mathrm{\φ}\left(k\right)=[\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{li}}^T{H}_{\mathrm{lc}}}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}{U}_c\left(k\right),{\mathrm{\λ}}_i\left(k\right)\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{li}}^T{H}_{\mathrm{lc}}}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}E]---\left(23\right)$

其中θ(k)为参数向量，φ(k)为回归向量，n(k)代表测量噪声的影响。由方程(21)用最小二乘法可估计出参数向量

，进而得出故障征兆参数。若考虑参数向量的动态模型为随机游走过程：

θ(k+1)＝θ(k)+ε(k)(24)

其中ε(k)为独立高斯随机向量，其均值为零，协方差矩阵为Q(k)，则可由(24)和(21)式用卡尔曼滤波方法来估计参数向量。计算公式如下：

K(k)＝P(k-1)φ^T(k)[R(k)+φ(k)P(k-1)φ^T(k)]^-1    (25)

$\widehat{\mathrm{\θ}}\left(k\right)=\widehat{\mathrm{\θ}}(k-1)+K\left(k\right)[r_i\left(k\right)-\mathrm{\φ}\left(k\right)\widehat{\mathrm{\θ}}(k-1)]---\left(26\right)$

P(k)＝P(k-1)+Q(k-1)-K(k)φ(k)P(k-1)-K(k)φ(k)Q(k-1)(27)

其中K(k)、P(k)为增益阵和协方差阵。

传感器残差计算方法：

假设非线性系统中执行器工作正常，仅考虑传感器故障。

对式(8)所表示的时间冗余测量方程，其在k时刻第l个工作点处的全解耦奇偶方程如下：

$r^l\left(k\right){=v}_l^T[Z\left(k\right)-H_{\mathrm{lc}}{U}_c\left(k\right)]---\left(28\right)$

其中，v_l为满足定义3的奇偶向量，r^l(k)为包含故障征兆信息的残差。Z(k)为传感器输出，当传感器正常时有Z(k)＝Y(k)。为了使残差r^l(k)仅对特定故障敏感，式中v_l应满足：

$V^{*}=\left\{v_l\right|v_l^T{H}^{*}=0\}---\left(29\right)$

其中H^*为全解耦奇偶空间矩阵，可用来使残差仅对特定传感器敏感，而对系统状态、扰动及其它传感器解耦。因此，构造矩阵H^*是获得对特定传感器敏感的全解耦奇偶方程的关键。

对于正常系统，r^l(k)可从式(8)得到

$r^l\left(k\right)=v_l^T[Y\left(k\right)-H_{\mathrm{lc}}{U}_c\left(k\right)]$

$=v_l^T[H_{l0}x(k-s)+H_{\mathrm{lc}}U\left(k\right)+H_{\mathrm{lw}}W\left(k\right)-H_{\mathrm{lc}}{U}_c\left(k\right)]---\left(30\right)$

$=v_l^T[H_{l0}x(k-s)+H_{\mathrm{lw}}W\left(k\right)]$

由上式可见，若v_l满足

则系统正常时r^l(k)为零。

为了使式(28)产生的残差仅对特定传感器敏感，可考虑残差对一个传感器解耦，而对其它传感器敏感，由式(8)可得

Y_j(k)＝H_jl0x(k-s)+H_jlcU(k)+H_jlwW(k)(31)

其中Y_j(k)＝[y(k-s)…y(k)]^T，为令对应传感器j输出为零时的系统的输出，矩阵H_jl0、H_jlc、H_jlw分别为将H_l0、H_lc、H_lw中的C_l、D_l、G_l用与传感器j对应的C_jl、D_jl、G_jl(j＝1，2，…，q)代替后形成。C_jl、D_jl、G_jl为将C_l、D_l、G_l中的第j行用0代替后得到的矩阵，例如，对应第1个传感器有

其中c_jl、d_jl、g_jl为矩阵C_jl、D_jl、G_jl中的行向量。当系统工作在工作点l处时，根据式(31)获得的残差将对传感器j不敏感。

由式(28)可知，在第l个工作点对传感器j不敏感的残差可表示为：

$r_j^l\left(k\right)=v_{\mathrm{jl}}^T[Z_j\left(k\right)-H_{\mathrm{jlc}}{U}_c\left(k\right)]---\left(32\right)$

Z_j(k)为对应Y_j(k)的实际输出。奇偶向量v_jl须与状态和扰动输出解耦，即满足：

$v_{\mathrm{jl}}^T=\left[\begin{array}{cc}{H}_{\mathrm{jl}0}& H_{\mathrm{jlw}}\end{array}\right]=0---\left(33\right)$

则系统对传感器j不敏感的残差为：

$r_j\left(k\right)=\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{jl}}^T[Z_j\left(k\right)-H_{\mathrm{jlc}}{U}_c\left(k\right)]}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}---\left(34\right)$

假设矩阵[H_jl0H_jlw]中不相关的列向量数为n_x，则式(33)有非零解(即奇偶向量存在)的充要条件为：

$s>\frac{n_x}{q}-1---\left(35\right)$

其中，s为数据窗的宽度。

传感器故障征兆估计方法：

传感器j的故障模型可表示为：

z_j(k)＝y^j(k)+f_j(k)(36)

其中y^j为传感器j的正常输出，z_j(k)为其实际输出。f_j(k)为传感器故障。传感器正常时f_j(k)＝0。表示为矩阵形式：

z(k)＝y(k)+f(k)(37)

其中：f(k)＝[f₁(k)f₂(k)…f_q(k)]^T。

假设在数据窗s+1内故障模型不变，则有：

Z(k)＝Y(k)+I^*f(k)(38)

其中Z(k)＝[z^T(k-s)z^T(k-s+1)…z^T(k)]^T为实际传感器输出，I^*＝[I₀I₀…I₀]^T为(s+1)q×q维矩阵，I₀是q×q维单位阵。

由(38)式可知：

Z_j(k)＝Y_j(k)+I^*f^(j)(k)(39)

其中Y_j(k)和f^(j)(k)分别为传感器j的输出和故障f_j用0代替后得到的系统输出和故障向量。

将式(39)代入式(34)，可得对j个传感器不敏感的残差为

$r_j\left(k\right)=\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{jl}}^T\left[(H_{\mathrm{jl}0}x(k-s)+H_{\mathrm{jlc}}U\left(k\right)+H_{\mathrm{jlw}}W\left(k\right)+I^{*}{f}^{\left(j\right)}\left(k\right)-H_{\mathrm{jlc}}{U}_c\left(k\right))\right]}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}---\left(40\right)$

由于

U(k)＝U_c(k)，所以：

$r_j\left(k\right)=\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_{\mathrm{jl}}^T{I}^{*}{f}^{\left(j\right)}\left(k\right)}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}=-{\mathrm{\φ}}_j\left(k\right){\mathrm{\θ}}_j\left(k\right)--\left(41\right)$

针对每一个传感器设计一个方程，可得

$r\left(k\right)=\frac{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)\·v_l^T{I}^{*}f\left(k\right)}{\underset{l=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(\mathrm{\ψ}\left(k\right)\right)}=\mathrm{\φ}\left(k\right)f\left(k\right)---\left(42\right)$

将模型误差等因素造成的不精确考虑为残差的噪声项，则有：

r(k)＝φ(k)f(k)+n(k)

其中r(k)＝[r₁(k)r₂(k)…r_q(k)]^T，v_l＝[v_1lv_2l…v_ql]^T，

f(k)＝[f₁(k)f₂(k)…f_q(k)]^T。

则f(k)为故障征兆向量，φ(k)为回归向量，n(k)是均值为零、协方差为R(k)的白噪声。用递推最小二乘方法或卡尔曼滤波方法得到故障征兆向量

步骤S2，采用AR模型对S1中产生的偏差序列进行建模，给出偏差预测值。具体实现方法如下：

AR(p)模型可表示为

式中

为常数，ε_t为纯随机过程，x(t-1)为时间序列t-1时刻的输出值，……，x(t-p)为序列t-p时刻的输出值。

由于AR模型的参数估计采用的是固定的一批数据，不能反映系统数据更新的特点，故采用滚动数据窗的方法给残差建模。随着数据窗的不断更新，模型的最优阶数也随之变化。确定模型的标准是使线性预报的方差达到最小，可采用FPE(Final Prediction Error)准则进行模型定阶。模型参数

可采用最小二乘法求取。模型建立成功后其预报即是线性最小方差意义下的平稳预报。

模型建立后，还要对模型进行平稳性检验，对模型残差进行白噪声检验，以确定建模是否成功。

在实际预测中，时间序列{x(t_l-N+1)，x(t_l-N+2)，…，x(t_l)}有可能是不平稳的，但因故障征兆一股变化较为缓慢，因此此不平稳的残差序列经过一阶差分有可能是平稳的。一阶差分公式为

$\left\{\begin{array}{c}\mathrm{\Δx}\left(0\right)=x\left(0\right)\\ \mathrm{\Δx}\left(t\right)=x\left(t\right)-x(t-1)\end{array}\right.---\left(44\right)$

当差分序列{Δx(t_l-N+1)，Δx(t_l-N+2)，…，Δx(t_l)}平稳时，即可利用AR模型对残差差分序列进行建模。

设t时刻对t+k时刻Δx(t+k)的估计值为

由式(44)可得，

$\left\{\begin{array}{c}\mathrm{\Δ}\hat{x}(t+1/t)=\hat{x}(t+1/t)-\hat{x}\left(t\right)\\ \hat{x}(t+k/t)=\hat{x}(t+k/t)-\hat{x}(t+k-1/t),(k\≥2)\end{array}\right.---\left(45\right)$

对t+k时刻进行故障预测，就是判断t+k时刻x(t+k)是否超出给定阈值。因此需要根据

进一步得到

将第1步到第k步的预测结果进行累加

$\mathrm{\Δ}\hat{x}(t+1/t)+\mathrm{\Δ}\hat{x}(t+2/t)+...+\mathrm{\Δ}\hat{x}(t+k/t)$

$=\hat{x}(t+1/t)-{\mathrm{\θ}}_e\left(t\right)+...+\hat{x}(t+k/t)-\hat{x}(t+k-1/t)$

$\hat{x}(t+k/t)=x\left(t\right)+\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\mathrm{\Δ}\hat{x}(t+i/t)---\left(46\right)$

即为在t时刻预测t+k时刻的残差值。

步骤S3，由偏差预测值结合其统计规律计算执行器或传感器故障发生概率，并用预测置信因子对预测准确性进行评价。具体实现方法如下：

故障发生概率计算方法：

设AR模型对的k步预测误差为ε_t(k)

${\mathrm{\ϵ}}_t\left(k\right)=\mathrm{\Δx}(t+k)-\mathrm{\Δ}\hat{x}(t+k/t)---\left(47\right)$

t时刻x(t+k)的预测误差为e_t(k)

$e_t\left(k\right)=x(t+k)-\hat{x}(t+k/t)---\left(48\right)$

可以证明，对x(t+k)的k步预测误差e_t(k)等于差分序列{Δx(t-N+1)，Δx(t-N+2)，…，Δx(t)}的1步预测到k步预测误差之和，即：

$e_t\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\ϵ}}_t\left(i\right)---\left(49\right)$

定理2：对于线性最小方差预报，其k步预报方差σ_t ²(k)可表示为

${\mathrm{\σ}}_t^2\left(k\right)={\mathrm{\σ}}_{\mathrm{\ϵ}}^2(1+G_1^2+...+G_{k-1}^2)---\left(50\right)$

为式(43)中模型参数

的估计值，R(0)，…，R(p)是数据序列的样本协方差函数在不同迟后时的值。G₁，…，G_k-1为格林函数

$G_i=-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{a}_j{G}_{i-j},$ G₀≡1，

R(j)＝E[(x(t)-μ)(x(t+j)-μ)]

其中μ为待预测序列的均值。当ε_t(k)服从均值为0方差为σ_t ²(k)的正态分布，e_t(k)也服从均值为0的正态分布，其方差

${\mathrm{\σ}}_{\mathrm{et}}^2\left(k\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\σ}}_t^2\left(i\right)---\left(51\right)$

设t+k时刻故障发生概率为α_t(k)，则

α_t(k)＝P{x(t+k)＞T_D}(52)

将式(48)代入式(52)可得

${\mathrm{\α}}_t\left(k\right)=P\{e_t\left(k\right)>T_D-\hat{x}(t+k)\}---\left(53\right)$

则

${\mathrm{\α}}_t\left(k\right)=\mathrm{\Φ}\left(\frac{\hat{x}(t+k)-T_D}{{\mathrm{\σ}}_{\mathrm{et}}\left(k\right)}\right)---\left(54\right)$

上式中

为t时刻对t+k时刻的预测值，T_D为诊断阈值，σ_et(k)为k步的预测误差。这也就是说，当预测结果恰好等于阈值时，预测的故障发生概率为50％；当预测的结果远大于阈值时，预测的故障发生概率为1；预测值很小时，预测的故障发生概率趋近于0。

预测置信因子计算方法：

定义1：系统故障：当t+k时刻执行器偏差超出阈值时，即认为此时卫星姿态控制系统发生故障，用C_t+k表示

Ct_+k＝{x(t+k)：x(t+k)≥T_D}(55)

其中x(t+k)为t+k时刻残差值，T_D为阈值。

定义2：预测系统故障：当t时刻预测系统在t+k时刻的残差超出阈值时，就预测k步后系统将发生故障，用A_t，k表示，则

$A_{t,k}=\{\hat{x}(t+k/t):\hat{x}(t+k/t)\≥T_D\}---\left(56\right)$

其中

为t时刻对t+k时刻残差的预测值。

定义3：预测正确率(probability of correct alarms)是事件A_t，k发生(即t时刻预测t+k时刻发生故障)的条件下事件C_t+k发生(即t+k时刻发生故障)的概率，即P{C_t+k|A_t，k}。

定义4：故障正确预测率(probability of detecting a fault)是在事件C_t+k已经发生的情况下，倒推t时刻事件A_t，k发生的概率，即P{A_t，k|C_t+k}。

定义5：无误预测率(probability of detecting no fault)是在事件C_t+k ^*(事件C_t+k的补集)已经发生的情况下，倒推t时刻事件A_t，k ^*(事件A_t，k的补集)发生的概率，即P{A_t，k ^*|C_t+k ^*}。

预测正确率与故障正确预测率是故障预测的两个主要评价参数。当事件A_t，k发生时，预测正确率P{C_t+k|A_t，k}与式(54)中的α_t(k)相等。故障正确预测率P{A_t，k|C_t+k}是一个后验概率，其大小与所建模型的精确程度及预测步长均有关系。由于残差序列本身具有不确定性，使得对建模误差的估计较为困难。当通过平稳性检验和模型残差的白噪声检验后，则认为建模成功，从而可忽略建模误差的影响。本发明主要考虑预测步长对故障预测准确性的影响。

定义6：预测置信因子(confidence factor)：是用来评价故障预测准确度的指标，反映由于预测步长增加等因素影响而导致的预测准确度的降低程度。k步预测的置信因子记为c_k。

显然，c_k与P{A_t，k|C_t+k}相关。根据残差x(t+k)的不同，分两种情况进行讨论。

1)x(t+k)≥T_D

此时事件C_t+k发生，

$P\left\{A_{t,k}\right|C_{t+k}\}$ (57)

$=P\{\hat{x}(t+k)>T_D|x(t+k)>T_D\}$

将式(48)代入式(57)，可得

P{A_t，k|C_t+k}

＝P{x(t+k)-e_t(k)＞T_D|x(t+k)＞T_D}(58)

＝P{e_t(k)＜x(t+k)-T_D|x(t+k)＞T_D}

设

ξ_e(t+k)＝|x(t+k)-T_D|(59)

则

P{A_t，k|C_t+k}＝P{e_t(k)＜ξ_e(t+k)}(60)

即当t+k时刻系统故障时，t时刻的故障正确预测率可转化为残差的k步预测误差e_t(k)小于ξ_e(t+k)的概率。

2)x(t+k)＜T_D

此时

$P\left\{A_{t,k}^{*}\right|C_{t+k}^{*}\}$

$=P\{x(t+k)-e_t\left(k\right)<T_D|x(t+k)<T_D\}---\left(61\right)$

$=P\{-e_t\left(k\right)<-(x(t+k)-T_D)|x(t+k)<T_D\}$

ξ_e(t+k)＝|x(t+k)-T_D|＝-(Δx(t+k)-T_D)

$P\left\{A_{t,k}^{*}\right|C_{t+k}^{*}\}=P\{-e_t\left(k\right)<{\mathrm{\&xi;}}_e(t+k)\}---\left(62\right)$

$=P\{e_t\left(k\right)<{\mathrm{\&xi;}}_e(t+k)\}$

即当t+k时刻系统正常时，其在t时刻的无误预测率可转化为残差的k步预测误差e_t(k)小于ξ_e(t+k)的概率。

由式(60)、(62)可知，在事件C_t+k发生的条件下A_t，k发生的概率P{A_t，k|C_t+k}，与事件C_t+k ^*发生的条件下A_t，k ^*发生的概率P{A_t，k ^*|C_t+k ^*}可统一为

P{e_t(k)＜ξ_e(t+k)}(63)

当e_t(k)服从正态分布时

$P\{e_t\left(k\right)<{\mathrm{\&xi;}}_e(t+k)\}=\mathrm{\&Phi;}\left(\frac{{\mathrm{\&xi;}}_e(t+k)}{{\mathrm{\&sigma;}}_{\mathrm{et}}\left(k\right)}\right)---\left(64\right)$

令

$P_r\left(k\right)=\mathrm{\&Phi;}\left(\frac{{\mathrm{\&xi;}}_e(t+k)}{{\mathrm{\&sigma;}}_{\mathrm{et}}\left(k\right)}\right)---\left(65\right)$

则P_r反映了系统预测的可信程度。当t+k时刻系统故障时，P_r与故障正确预测率相等；当t+k时刻系统正常时，P_r与无误预测率相等。说明P_r越大，则系统正确预测的可能性越大。因此可令预测置信因子为

c_k＝P_r(k)(66)

注意到σ_et(k)随着k的增大而增大，且

由一股经验可知，当k→∞时，其置信度应为0。因此对c_k稍作调整，令

c_k＝P_r(k)-0.5   (67)

此时c_k的范围为(0，0.5)。通常我们习惯于置信因子在(0，1)范围内。将c_k归一化，得

$c_k=\frac{P_r\left(k\right)-0.5}{0.5}=2\mathrm{\&Phi;}\left(\frac{{\mathrm{\&xi;}}_e(t+k)}{{\mathrm{\&sigma;}}_{\mathrm{et}}\left(k\right)}\right)-1---\left(68\right)$

上式说明c_k与e_t(k)和ξ_e(t+k)有关。而ξ_e(t+k)在t时刻是未知的，确定为ξ_e(t+k)＝3σ_et(1)。

本发明基于现有故障预测方法的不足，针对工程实际中普遍存在的非线性系统，给出了基于模糊奇偶方程和模糊基函数网络的故障预测方法，并针对预测的不准确性，以故障发生概率的形式给出了预测结果，能够给人以更直观的印象，并以预测置信因子来反映了预测结果的准确程度。

以上的实施例仅是对本发明的优选实施方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通工程技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明的权利要求书确定的保护范围内。

Claims (4)

1.一种基于模糊奇偶方程及AR模型的非线性系统故障预测方法，其特征在于，包括以下步骤：

采用模糊奇偶方程方法估计非线性系统执行器或传感器偏差；

采用AR模型对所述执行器、传感器产生的偏差序列进行建模，给出偏差预测值；

由偏差预测值结合其统计规律计算执行器或传感器故障发生概率，并用预测置信因子对预测准确性进行评价，其中：

所述采用模糊奇偶方程方法估计非线性系统执行器、传感器偏差的步骤还包括：

对非线性系统在各工作点处线性化，得到局部线性模型；

离线计算出各局部线性模型的奇偶方程；

利用模糊奇偶方程方法计算出系统的残差，根据残差进行故障诊断；

所述利用模糊奇偶方程方法计算系统的残差，根据残差进行故障诊断步骤还包括：执行器残差的计算和故障征兆评估，以及传感器残差的计算和故障征兆评估；

所述采用AR模型对所述执行器、传感器产生的偏差序列进行建模的步骤中，采用滚动数据窗方法给残差建模，并采用FPE准则进行模型定阶；

所述预测置信因子用c_k表示，由残差的统计规律确定预测置信因子，当残差正态分布时，预测置信因子表示为

其中ξ_e(t+k)＝3σ_et(1)，σ_et(k)为k步的预测误差；

所述ξ_e(t+k)设为ξ_e(t+k)＝|x(t+k)-T_D|，其中x(t+k)为t+k时刻残差值，T_D为阈值。

2.根据权利要求1所述的基于模糊奇偶方程及AR模型的非线性系统故障预测方法，其特征在于，所述采用AR模型对所述执行器、传感器产生的偏差序列进行建模的步骤中还包括对模型进行平稳性检验，对模型残差进行白噪声检验的步骤。

3.根据权利要求1所述的基于模糊奇偶方程及AR模型的非线性系统故障预测方法，其特征在于，所述预测置信因子用于反映由于预测步长增加因素而导致的预测准确度的降低程度。

4.根据权利要求3所述的基于模糊奇偶方程及AR模型的非线性系统故障预测方法，其特征在于，所述故障预测包括预测正确率和故障正确预测率两个主要评价参数。