CN103761450A - 一种基于模糊自适应预测的动态过程故障预报方法 - Google Patents

Abstract

本发明属于化工生产过程故障诊断与预报技术领域，涉及一种基于模糊自适应预测的动态过程故障预报方法，先基于模糊自适应无迹卡尔曼预测方法依次确定故障参数和状态变量，再根据故障参数的标准差判断是否进行未来多步预测；然后计算未来多步预测值；最后利用测量输出变量的西格玛点计算故障预测概率，再根据故障预测概率进行故障预测；其整体工艺过程简单，原理可靠，计算参数准确，估计故障状态精度高，适应范围广，逻辑性强，环境友好。

Description

一种基于模糊自适应预测的动态过程故障预报方法

技术领域：

本发明属于化工生产过程故障诊断与预报技术领域，涉及一种基于模糊自适应预测的动态过程故障预报方法。

背景技术：

随着化学工业生产规模的不断扩大和生产过程的日益复杂，一旦发生过程故障就会影响产品质量，甚至造成人员伤亡和生态危机，为保证生产过程的稳定性、减小故障损失，在自动化生产过程中迫切希望预知过程未来运行状态，在过程超出正常运行区域之前尽可能早地检测到故障，将故障消灭在萌芽状态。故障预报技术是利用故障发生初期的微小异常信息通过计算过程变量的未来预测值，对即将发生的异常事故提前报警，而且故障预报是设备状态维修的关键技术，有助于解决维修不足和检修过剩问题。因此，故障预报技术是提高过程安全、减小过程故障损失的关键。故障预报是故障检测和诊断技术的新兴研究领域；目前，相关研究成果较少，2001年Juricek和Seborg等提出预测监控概念，即在检测到故障之后利用过程变量预测值判断过程将来是否超出紧急控制限。目前，故障预报技术主要围绕两个方面开展研究：一是如何准确跟踪故障过程、准确预测未来过程信息；二是如何利用有限的未来预测信息判断过程的未来运行状态，其中，对于如何准确跟踪故障过程、准确预测未来过程信息的问题，亟需研究具有较强跟踪能力的过程状态和输出预测方法，现有的预测方法包括卡尔曼预测器、粒子预测器、强跟踪粒子滤波器、模糊自适应无迹（unscented）卡尔曼预测器、支持向量机和基于相关向量机模糊模型的预测方法；对于如何利用有限的未来预测信息判断过程的未来运行状态的问题，Juricek和Seborg等直接将未来输出预测值区间与控制限进行比较；胡昌华等利用强跟踪粒子滤波器估计过程状态变量，根据状态变量的粒子计算故障预测概率；模糊自适应unscented卡尔曼预测方法在论文（Tian Xuemin,Cao Yuping,ChenSheng,Process Fault Prognosis Using a Fuzzy-Adaptive Unscented Kalman Predictor,International Journal of Adaptive Control and Signal Processing,2011，25（9）：813–830/田学民，曹玉苹，陈生，基于模糊自适应无迹卡尔曼预测的过程故障预报，国际自适应控制与信号处理杂志，2011，25（9）：813–830,）中公开发表；但是，就上述现有技术而言，对于有的过程模型，状态变量没有实际物理意义，其正常运行区间难以确定，而且对于非线性动态过程，尚未见有逻辑完整、实用性较强的实时故障预报系统。

发明内容：

本发明的目的在于克服现有技术存在的缺点，寻求设计提供一种基于模糊自适应预测的动态过程故障预报方法，通过建立适用于非线性动态过程的实时故障预报系统，先基于模糊自适应unscented卡尔曼预测方法依次估计故障参数和状态变量，再根据故障参数的标准差判断是否进行未来多步预测；然后使用unscented卡尔曼预测方法计算未来多步预测值；最后利用测量输出变量的sigma点计算故障预测概率，再根据故障预测概率进行故障预测。

为了实现上述目的，本发明实现动态过程故障预报的具体过程包括以下步骤：

（1）、先采集故障参数的实时测量数据，使用模糊自适应无迹卡尔曼预测方法计算故障参数后再确定过程状态；

（2）、根据步骤（1）计算的故障参数的标准差判断是否进行未来多步预测；若进行未来多步预测，则进行步骤（3）；否则，返回执行步骤（1）；

（3）、使用模糊自适应无迹卡尔曼预测方法计算未来多步预测值，利用未来测量输出变量的西格玛（sigma）点计算故障预测概率；

（4）、根据故障预测概率判断动态过程将来是否超出正常运行区域；若超出正常运行区域，则进行故障预报警后返回执行步骤（1）；否则直接返回执行步骤（1）。

本发明涉及的模糊自适应无迹卡尔曼预测方法采用非线性动态系统模型，非线性动态系统模型为：

x(k+1)=f(x(k),u(k))+d(k)+w(k)   （1）

y(k)=H(k)x(k)+v(k)         （2）

其中k表示采样时刻，k为自然数，1≤k≤1000，x表示n维的状态向量，n为自然数，1≤n≤10，u表示l维的测量输入向量，l为自然数，1≤l≤10，d表示n维的不可测故障项，f(·)表示n维非线性函数向量，对x连续可导，y表示m维的测量输出向量，m为自然数，1≤m≤10，H表示m×n维的矩阵，w表示n维的过程噪声，服从均值为零协方差矩阵为Q的正态分布，Q表示n×n维的协方差矩阵，v表示m维的测量噪声，服从均值为零协方差矩阵为R的正态分布，R表示m×m维的协方差矩阵，w与v不相关；u(k)∈R^l×1，x(k)∈R^n×1，y(k)∈R^m×1分别表示测量输入、状态和测量输出；d(k)∈R^n×1表示不可测故障项；w(k)∈R^n×1，v(k)∈R^m×1分别表示过程噪声和测量噪声，二者不相关且分别服从均值为零、方差为Q(k)和R(k)的正态分布；当过程正常运行时，故障项d(k)=0；当过程出现故障时，故障项d(k)≠0；故障预报针对缓慢变化型故障采用的故障模型为：

d(k)=Γ(k,k₀)θ       （3）

其中θ表示2n维的故障参数向量，k₀表示故障起始时间，Γ(k,k₀)为n×2n维的时间函数矩阵；以2阶系统为例，故障模型为：

$d\left(k\right)=\mathrm{\Γ}(k,k_0)\mathrm{\θ}=\left[\begin{array}{c}{\mathrm{\θ}}_{\mathrm{1,0}}+{\mathrm{\θ}}_{1,s}{T}_s(k-k_0)\\ {\mathrm{\θ}}_{\mathrm{2,0}}+{\mathrm{\θ}}_{2,s}{T}_s(k-k_0)\end{array}\right]---\left(4\right)$

其中

$\mathrm{\Γ}(k,k_0)=\left[\begin{array}{cccc}1& 0& T_s(k-k_0)& 0\\ 0& 1& 0& T_s(k-k_0)\end{array}\right]$

θ=[θ_1,0θ_2,0θ_1,sθ_2,s]^T

其中T_s为采样周期，0≤T_s≤1分钟，θ_h,0表示故障参数向量θ的第h个元素，θ_h,s表示故障参数向量θ的第2+h个元素，h=1,2；

故障参数θ的动态模型为

$\left\{\begin{array}{c}\mathrm{\θ}(k+1)=\mathrm{\θ}\left(k\right)+\mathrm{\ϵ}\left(k\right)\\ y(k+1)=f(x\left(k\right),u\left(k\right))+\mathrm{\Γ}(k,k_0)\mathrm{\θ}(k+1)+w^{\′}\left(k\right)\end{array}\right.---\left(5\right)$

其中ε表示2n维的过程噪声，服从均值为零协方差矩阵为Q_θ的正态分布，Q_θ表示2n×2n维的协方差矩阵，w′(k)=w(k)-Γ(k,k₀)ε(k)，w′表示m维的测量噪声，服从均值为零协方差矩阵为R_θ的正态分布，R_θ表示m×m维的协方差矩阵。

本发明涉及的模糊自适应无迹卡尔曼预测方法的工艺步骤包括单步预测、模糊自适应更新和未来多步预测三个步骤，其具体过程如下：

（A1）、单步预测：对公式（1）（2）所示的非线性动态系统模型，构造增广状态向量x_a(k-1)=[x^T(k-1)  w^T(k-1)  v^T(k-1)]^T，设x_a(k-1)的第i个西格玛（sigma）点为X_i(k-1)，x_a(k-1)的第i个西格玛（sigma）点的权值为ω_i，X_i(k-1)为n_a维向量，0≤i≤2n_a，n_a=2n+m，a表示增广，X_i(k-1)和ω_i的计算公式为：

$X_i(k-1)=\left\{\begin{array}{cc}{\tilde{x}}_a(k-1),& i=0\\ {\tilde{x}}_a(k-1)+A_i(k-1),& i=1,...,{2n}_a\\ {\tilde{x}}_a(k-1)-A_{i-n_a}(k-1),& i=n_a+1,...,{2n}_a\end{array}\right.{\mathrm{\ω}}_i=\left\{\begin{array}{cc}\frac{\mathrm{\η}}{n_a+\mathrm{\η}},& i=0\\ \frac{1}{2(n_a+\mathrm{\η})},& i=1,...,{2n}_a\end{array}\right.---\left(6\right)$

其中 ${\tilde{x}}_a(k-1)={\left[\begin{array}{ccc}{\hat{x}}^T(k-1)& 0_{1\×n}& 0_{1\×m}\end{array}\right]}^T,$

表示第k-1个采样时刻状态x的估计值，0_1×n表示1行n列的零矩阵，0_1×m表示1行m列的零矩阵，向量A_i表示矩阵A的第i列，矩阵A为n_a×n_a维，满足A(k-1)A^T(k-1)=(n_a+η)P_a(k-1)，矩阵P_a(k-1)为n_a×n_a维，表示第k-1个采样时刻增广状态的误差协方差矩阵，参数η根据状态x的分布确定，对于高斯分布η=3-n_a；记预测状态的第i个西格玛（sigma）点为

状态预测值为

其误差协方差矩阵为P(k|k-1)，其中为n维向量，

表示在第k-1个采样时刻预测的第k个采样时刻的x的估计值，P(k|k-1)为n×n维矩阵，表示在第k-1个采样时刻预测的第k个采样时刻的误差协方差矩阵，则

和P(k|k-1)的计算公式如下：

${\tilde{X}}_i\left(k\right)=f({X_i}^x(k-1),u(k-1))+{X_i}^w(k-1)+\mathrm{\Γ}(k-1,k_0)\widehat{\mathrm{\θ}}(k-1)---\left(7\right)$

$\hat{x}\left(k\right|k-1)=\underset{i=0}{\overset{{2n}_a}{\mathrm{\Σ}}}{\mathrm{\ω}}_i{\tilde{X}}_i\left(k\right)---\left(8\right)$

$P\left(k\right|k-1)=\underset{i=0}{\overset{{2n}_a}{\mathrm{\Σ}}}{\mathrm{\ω}}_i[{\tilde{X}}_i\left(k\right)-\hat{x}\left(k\right|k-1)]{[{\tilde{X}}_i\left(k\right)-\hat{x}\left(k\right|k-1)]}^T---\left(9\right)$

其中

为故障参数θ(k-1)的估计值，X_i ^x(k-1)和X_i ^w(k-1)分别表示X_i(k-1)中与状态和过程噪声对应的部分，X_i ^x(k-1)由X_i(k-1)的第1至n个元素组成，X_i ^w(k-1)由X_i(k-1)的第n+1至2n个元素组成；记预测输出的第i个西格玛（sigma）点为

输出预测值为

其中

为m维向量，

表示在第k-1个采样时刻预测的第k个采样时刻的y的估计值，则和

的计算公式为：

${\tilde{Y}}_i\left(k\right|k-1)=H\left(k\right){\tilde{X}}_i\left(k\right)+{X_i}^v(k-1)---\left(10\right)$

$\hat{y}\left(k\right|k-1)=\underset{i=0}{\overset{{2n}_a}{\mathrm{\Σ}}}{\mathrm{\ω}}_i{\tilde{Y}}_i\left(k\right|k-1)---\left(11\right)$

其中X_i ^v(k-1)表示X_i(k-1)中与测量噪声对应的部分，X_i ^v(k-1)由X_i(k-1)的第2n+1至2n+m个元素组成；

（A2）、模糊自适应更新：利用输出预测误差及其协方差矩阵的估计值

计算自适应调节参数λ(k)，

$P_{\tilde{y}}\left(k\right)=\left\{\begin{array}{cc}V\left(k\right),& k=1\\ \frac{\mathrm{\ρV}(k-1)+V\left(k\right)}{1+\mathrm{\ρ}},& k\≥2\end{array}\right.---\left(12\right)$

$\mathrm{\λ}\left(k\right)=\max \{\frac{\mathrm{tr}[P_{\tilde{y}}\left(k\right)-\mathrm{\ζR}\left(k\right)]}{\mathrm{tr}\left[H\left(k\right)P\left(k\right|k-1)H^T\left(k\right)\right]},1\}---\left(13\right)$

其中

ρ为遗忘因子，0≤ρ≤1，ζ为弱化因子，1≤ζ≤10，tr[·]为求迹算子；由于过程故障运行状态的模糊性，利用模糊逻辑系统计算自适应调节参数λ^*(k)。令模糊集合“Normal”和“Failure”分别表示正常运行状态和故障运行状态，对应的隶属度函数分别为μ_N(y(k))和μ_F(y(k))，且μ_N(y(k))+μ_F(y(k))=1；模糊规则为：

R¹:IF y(k)is“Normal”,THEN

R²:IF y(k)is“Failure”,THEN

设测量输出向量y(k)=[y₁(k),…y_z(k),…ym(k)]^T服从均值为μ协方差为Σ的正态分布，其中y_z(k)为y(k)的第z个元素，1≤z≤m，μ为m维向量，Σ为m×m维矩阵，则

$\begin{array}{c}{\mathrm{\&mu;}}_F\left(y\left(k\right)\right)=P(-|y_1\left(k\right)|\&le;{\mathrm{\&psi;}}_1<|y_1\left(k\right)|,\&CenterDot;\&CenterDot;\&CenterDot;-|y_z\left(k\right)|\&le;{\mathrm{\&psi;}}_z<|y_z\left(k\right)|,\&CenterDot;\&CenterDot;\&CenterDot;,-|y_m\left(k\right)|\&le;{\mathrm{\&psi;}}_m<|y_m\left(k\right)\left|\right)\\ =\frac{1}{\sqrt{{\left(2\mathrm{\&pi;}\right)}^m\left|\mathrm{\&Sigma;}\right|}}{\&Integral;}_{-\left|y_1\left(k\right)\right|}^{\left|y_1\left(k\right)\right|}\&CenterDot;\&CenterDot;\&CenterDot;{\&Integral;}_{-\left|y_z\left(k\right)\right|}^{\left|y_z\left(k\right)\right|}\&CenterDot;\&CenterDot;\&CenterDot;{\&Integral;}_{-\left|y_m\left(k\right)\right|}^{\left|y_m\left(k\right)\right|}{e}^{-\frac{1}{2}{(\mathrm{\&psi;}-\mathrm{\&mu;})}^T{\mathrm{\&Sigma;}}^{-1}(\mathrm{\&psi;}-\mathrm{\&mu;})}\mathrm{d\&psi;}\end{array}---\left(14\right)$

其中向量ψ=[ψ₁,…ψ_z,…,ψ_m]^T，ψ_z为ψ的第z个元素，P(·)表示概率函数，利用Monte-Carlo方法或分区自适应方法等数值积分算法计算式（14）；先计算自适应调节参数λ^*(k)

${\mathrm{\&lambda;}}^{*}\left(k\right)={\mathrm{\&mu;}}_N\left(y\left(k\right)\right){\mathrm{\&lambda;}}_1^{*}\left(k\right)+{\mathrm{\&mu;}}_F\left(y\left(k\right)\right){\mathrm{\&lambda;}}_2^{*}\left(k\right)---\left(15\right)$

再利用λ^*(k)调节P(k|k-1)

$\stackrel{\&OverBar;}{P}\left(k\right|k-1)={\mathrm{\&lambda;}}^{*}\left(k\right)P\left(k\right|k-1)---\left(16\right)$

其中矩阵

表示矩阵P(k|k-1)的模糊自适应修正值，则卡尔曼增益K(k)、状态修正值及其误差协方差矩阵P(k)的计算公式为：

$K\left(k\right)=\stackrel{\&OverBar;}{P}\left(k\right|k-1)H^T\left(k\right){[H\left(k\right)\stackrel{\&OverBar;}{P}\left(k\right|k-1)H^T\left(k\right)+R\left(k\right)]}^{-1}---\left(17\right)$

$\hat{x}\left(k\right)=\hat{x}\left(k\right|k-1)+K\left(k\right)(y\left(k\right)-\hat{y}\left(k\right|k-1))---\left(18\right)$

$P\left(k\right)=\stackrel{\&OverBar;}{P}\left(k\right|k-1)-K\left(k\right)[H\left(k\right)\stackrel{\&OverBar;}{P}\left(k\right|k-1)H^T\left(k\right)+R\left(k\right)]K^T\left(k\right)---\left(19\right);$

（A3）、未来多步预测；在第k个采样时刻，重复迭代执行步骤（A1）得到未来p步的输出预测值

表示在第k个采样时刻预测的第k+j个采样时刻的y的预测值，j=1,2,…,p，p为未来预测步数，1≤p≤100，其中

P(k+j)=P(k+jk)（20）

$\hat{x}(k+j)=\hat{x}(k+j|k)---\left(21\right)$

${\tilde{Y}}_i(k+j|k)=H(k+j){\tilde{X}}_i(k+j)+{X_i}^v(k+j-1)---\left(22\right)$

P(k+j|k)表示在第k个采样时刻预测的第k+j个采样时刻的误差协方差矩阵，

表示在第k个采样时刻预测的第k+j个采样时刻的x的估计值，

表示在第k个采样时刻预测的第k+j个采样时刻的西格玛（sigma）点。

本发明涉及的步骤（1）确定过程状态时，在每个采样时刻，先计算故障参数的估计值

然后确定状态x，其具体过程为：

（1.1）、在第k个采样时刻，根据故障参数的动态模型公式（5），将故障参数作为状态，利用模糊自适应无迹卡尔曼预测方法的步骤（A1）和步骤（A2）计算故障参数的估计值

或利用扩展卡尔曼滤波器估计；

（1.2）、根据过程动态模型公式（1）和（2），利用模糊自适应无迹卡尔曼预测方法的步骤（A1）和步骤（A2）计算过程状态的估计值

本发明涉及的步骤（2）中根据故障参数的标准差判断是否进行未来多步预测，其判断准则为：在第k个采样时刻，若

且…且

且…且

其中std(·)表示标准差函数，σ_h表示标准差的阈值，根据过程噪声w的标准差确定，h=1,2,…,n，则进行未来预测的采样时刻k_p=k，执行步骤（3）计算未来p步预测值；否则，返回执行步骤（1），继续跟踪过程的变化趋势；其中阈值σ_h根据过程噪声w的标准差确定。

本发明涉及的步骤（3）中未来多步预测值和故障预测概率的计算方法为：

（3.1）、在第k_p个采样时刻，利用模糊自适应unscented卡尔曼预测方法的步骤（A3）计算未来p步的输出预测值

及其西格玛（sigma）

其中j=1,2,…,p，i=0,1,2,…,2n_a；

（3.2）、对未来输出预测值的西格玛（sigma）进行分析，若

则则第i个故障指标函数F_i(k_p+j|k_p)=1；否则F_i(k_p+j|k_p)=0；其中

为第i个西格玛（sigma）的第z个元素；Bz为第z个输出变量的故障区间，由输出变量的控制限决定；i=0,1,2,…,2n_a,z=1,2,…,m；按照公式（24）计算故障预测概率：

$P_f(k_p+j|k_p)=\underset{i=0}{\overset{{2n}_a}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_i{F}_i(k_p+j|k_p)---\left(24\right).$

本发明涉及的步骤（4）中故障预报警判断准则为：在第k_p个采样时刻，根据P_f(k_p+j|k_p)判断过程将来是否超出正常运行区域发生故障，若P_f(k_p+q|k_p)=1且P_f(k_p+q+g|k)=1，则预报警采样时刻k_b=k_p，发出故障预报警，然后返回执行步骤（1）；否则，直接返回执行步骤（1），继续跟踪过程变化趋势；其中q=1,2,…,p-1,g=1,2,…,p-q。

本发明利用模糊自适应无迹卡尔曼预测方法进行故障预测的原理是：故障发生时，过程由平稳运行状态突变到故障状态，对滤波器的跟踪能力要求较高；结合故障状态的模糊性和强跟踪滤波，利用模糊逻辑系统计算自适应调节参数，在线调节状态估计误差协方差矩阵，强迫滤波器保持对实际系统的跟踪；过程处于故障运行状态时，输出预测误差是由故障模型失配引起的，先计算故障参数后估计过程状态，能够提高未来预测的准确性。

本发明利用测量输出变量的sigma点计算故障预测概率的原理是：对于有的非线性动态过程模型，状态变量没有实际物理意义，其正常运行区间难以确定，动态过程的测量输出变量具有实际物理意义，而且操作员根据生产运行状态对测量输出变量设置控制限，根据控制限能够确定变量的正常区间和故障区间。

本发明设置未来预测判断准则的原理是：每个采样时刻均进行未来多步预测费时费力且占用较大存储空间，在过程处于正常运行状态时进行故障预报没有实际意义；当过程处于故障运行状态时，若故障参数估计值不收敛，未来多步预测值不可信，进行故障预报会导致错误的报警，因此，要判断过程是否需要进行未来多步预测，通过分析故障参数估计值

发现：故障状态下

的标准差明显小于正常状态下的标准差，利用

序列的标准差信息确定是否进行未来多步预测。

本发明与现有技术相比具有以下优点：一是结合模糊逻辑系统和强跟踪滤波的优点，利用模糊自适应调节参数在线调节状态估计误差协方差矩阵，使得滤波器对故障过程的跟踪能力增强，先计算故障参数后确定过程状态变量，有效提高未来预测信息的准确性，实现动态过程故障预报，适应非线性动态过程实时故障预报的需求；二是利用测量输出变量的sigma点计算故障预测概率，测量输出变量具有实际物理意义和操作控制限，可根据控制限确定变量的正常区间和故障区间，解决基于状态变量的方法缺乏实际物理意义和控制限的问题；三是设置未来预测判断准则，在过程发生故障且故障参数的估计值收敛时进行未来多步预测，有效减少不必要的预测计算量，降低盲目预测导致的误报率，该逻辑判断环节增强了实时故障预报系统的逻辑完整性和实用性；其整体工艺过程简单，原理可靠，计算参数准确，估计故障状态精度高，适应范围广，逻辑性强，环境友好。

附图说明：

图1为本发明的动态过程故障预报工作原理流程示意图。

图2为本发明实施例发生故障1时故障参数的估计结果比较曲线图，实线是基于本发明先估计故障参数后估计过程状态得到的估计结果，虚线是过程状态和故障参数联合估计得到的估计结果；其中（a）为θ_1,0的估计结果，（b）为θ_2,0的估计结果，（c）为θ_1,s的估计结果，（d）为θ_2,s的估计结果。

图3为本发明实施例发生故障1时输出变量未来预测值和实测值的比较曲线图；其中（a）为k_p=123时浓度C_A预测值和实测值的比较结果，（b）为k_p=124时浓度C_A预测值和实测值的比较结果，（c）为k_p=125时浓度C_A预测值和实测值的比较结果，（d）为k_p=123时温度T预测值和实测值的比较结果，（e）为k_p=124时温度T预测值和实测值的比较结果，（f）为k_p=125时温度T预测值和实测值的比较结果。

图4为本发明实施例发生故障1时的未来预测误差比较曲线图；其中（a）为k_p=123时的预测误差e²，（b）为k_p=124时的预测误差e²，（c）为k_p=125时的预测误差e²。

图5为本发明实施例发生故障1时的故障预测概率曲线图。

具体实施方式：

下面通过实施例并结合附图对本发明进行进一步说明。

实施例：

本实施例将动态过程故障预报方法应用到非绝热连续搅拌反应器（Continuous StirredTank Reactor,CSTR）过程中，非绝热连续搅拌反应器是一个典型的化工生产单元，反应器内发生一级不可逆放热反应A_r→B_r，热量通过冷却夹套带走，具体模型描述如下：

$\frac{{\mathrm{dC}}_A}{\mathrm{dt}}=-k_0{e}^{-\frac{E}{\mathrm{RT}}}{C}_A+\frac{F(C_F-C_A)}{V}$

$\frac{\mathrm{dT}}{\mathrm{dt}}=\frac{(-\mathrm{\&Delta;H})k_0{e}^{-\frac{E}{\mathrm{RT}}}{C}_A}{\mathrm{\&rho;}{C}_P}+\frac{F(T_F-T)}{V}+\frac{{\mathrm{UA}}_a(T_J-T)}{\mathrm{V\&rho;}{C}_P}$

其中A_a为传热面积，C_A为反应物A_r的浓度，C_F为进料浓度，C_p为反应物比热，E为反应活化能，F为进料流量，ΔH为单位反应热，k₀为反应频率因子，R为气体常数，t为时间，T为反应温度，T_F为进料温度，T_J为冷却剂出口温度，U为传热系数，V为反应器体积，ρ为反应物密度，反应器参数和常规操作条件如表1所示，非绝热连续搅拌反应器过程常发生的缓变型故障有进料温度变化、进料浓度变化、热交换器结垢等故障，其中，进料温度变化和进料浓度变化会影响输出产品的质量，热交换器结垢不仅影响反应温度和输出产品的质量，而且容易因为局部过热出现金属变形，影响热交换器的安全运行，模拟三种缓变型过程故障，具体故障描述如表2所示，基于Matlab软件平台对非绝热连续搅拌反应器正常运行过程和故障过程进行蒙特卡罗仿真，采样周期为0.2分钟，故障均从采样时刻k=101开始；

表1：CSTR参数和常规操作条件

表2：故障描述

定义过程状态和测量输出向量y=x=[C_A  T]^T，输入向量u=[T_J F]^T，则故障参数θ=[θ_1,0θ_2,0θ_1,sθ_2,s]^T；测量输出变量的控制限为C_A∈[0.05,0.15]，T∈[373,393]，所以输出变量故障区间B₁={C_A|C_A≤0.05或C_A≥0.15}，B₂={T|T≤373或T≥393}；模糊自适应预测方法参数设置：过程噪声参数Q(k)=diag([10^-7  10^-3])，Q_θ(k)=diag([10^-7  10-310^-7  10^-3])，测量噪声参数R(k)=R_θ(k)=diag([10^-6  0.04])，其中diag(·)表示对角阵创建函数；为控制调节速度，设置模糊自适应调节参数1≤λ^*(k)≤2；测量输出变量的概率分布参数μ=[0.1383.78]^T，

根据的数量级确定未来多步预测的判断阈值σ₁=10^-4,σ₂=10^-2；设置未来预测步数p=50；采用基于模糊自适应预测的动态过程故障预报方法对非绝热连续搅拌反应器进行实时故障预报的具体实施步骤为：

（1）、假设当前时刻为第k个采样时刻，采集实时测量数据u(k)，y(k)，先使用模糊自适应unscented卡尔曼预测方法估计故障参数，然后估计过程状态，具体流程为：

1.1）、根据故障参数的动态模型公式（5），将故障参数作为状态，利用模糊自适应unscented卡尔曼预测方法的公式（6）-（11）计算故障参数的预测值利用公式（12）-（19）计算故障参数的修正值

1.2）、根据过程的动态模型公式（1）和（2），利用模糊自适应unscented卡尔曼预测方法的公式（6）-（11）计算状态预测值

和输出预测值

利用公式（12）-（19）计算状态的修正值

（2）、根据故障参数的标准差判断是否进行未来多步预测，若

且

则进行未来预测的采样时刻k_p=k，执行步骤（3）计算未来50步的输出预测值。其中

为序列

的标准差，

的计算方法与此相同；否则，返回执行步骤（1），继续跟踪过程的变化趋势；

（3）、使用模糊自适应unscented卡尔曼预测方法计算未来多步预测值，利用未来测量输出变量的sigma点计算故障预测概率，具体流程为：

3.1）、对第kp个采样时刻，利用模糊自适应unscented卡尔曼预测方法的公式（6）-（9），（20）-（23）迭代计算未来50步的输出预测值

及其sigma点

其中j=1,2,…,50，i=0,1,2,…,12。

3.2）、对未来输出预测值的sigma点进行分析，若

或

则故障指标函数F_i(k_p+j|k_p)=1；否则F_i(k_p+j|k_p)=0；其中j=1,2,…,50，i=0,1,2,…,12；按照公式（24）计算故障预测概率函数P_f(k_p+j|k_p)；

（4）、根据故障预测概率判断过程将来是否超出正常运行区域发生故障，若P_f(k_p+q|k_p)=1且P_f(k_p+q+g|k)=1，则故障预报警采样时刻k_b=k_p，发出故障预报警，然后返回执行步骤（1）；否则，直接返回执行步骤（1），继续跟踪过程变化趋势；其中q=1,2,…,49,g=1,2,…,50-q。

本实施例为比较输出向量的预测效果，定义归一化平方预测误差e²及平均归一化平方预测误差（Averaged normalized square prediction errors,ASE），其中

$e^2(k_p+j)={\left|\right|{\mathrm{\&Sigma;}}^{-\frac{1}{2}}(y(k_p+j)-\hat{y}(k_p+j|k_p))\left|\right|}_2^2$

$\mathrm{ASE}=\frac{1}{50}\underset{j=1}{\overset{50}{\mathrm{\&Sigma;}}}{e}^2(k_p+j);$

本实施例涉及的故障1是进料温度T_F按照斜率2K/min缓慢上升，故障预报结果如图2-5所示，图2为故障参数θ=[θ_1,0θ_2,0θ_1,sθ_2,s]^T的估计结果，可以看出本发明方法得到的故障参数波动较小，收敛效果明显优于联合估计方法；利用步骤（2）判断得出在k_p=124时进行未来多步预测，表明故障参数收敛较快，在故障发生23步后开启未来预测，为对比预测结果，分别在k_p=123,124,125时进行未来50步预测，预测结果如图3所示，归一化平方预测误差e²如图4所示，平均归一化平方预测误差ASE如表3所示；从图3可以看出，与k_p=123时预测结果相比，k_p=124,125时的预测值更准确；从图4可以看出，k_p=123时的未来预测误差较大且随预测步长的增加而增大，这是由于k_p=123时故障参数未收敛，没必要进行未来多步预测；从表3可以看出，k_p=124,125时的ASE明显小于k_p=123时的计算结果；图5是在k_p=124时得到的故障预测概率，可见未来故障预测概率升至1，预测过程未来将超出正常运行区域，因此在k_b=124时发出故障预报警；由于温度的实际观测值在k_f=151时超出控制限，所以提前δk=k_f-k_b=27步预报；虽然在第124个采样时刻过程处于正常运行区域，预报警有利于操作员及时检查并采取维修措施将故障消灭在萌芽状态。

本实施例涉及的故障2是进料浓度C_F按照斜率0.01(mol/L)/min缓慢上升，利用步骤（2）判断得出在k_p=127时进行未来50步预测，k_p=126,127,128时的ASE如表3所示，由表3可见，k_p=127,128时的ASE明显小于k_p=126时的计算结果；利用步骤（4）判断得出在k_b=127时发出故障预警，由于温度的实际观测值在kf=154时超出控制限，所以提前δk=k_f-k_b=27步预报。

本实施例涉及的故障3是UA_a按照斜率125(J/(min·K))/min缓慢下降，模拟热交换器结垢故障，利用步骤（2）判断得出在k_p=130时进行未来50步预测，k_p=129,130,131时的ASE如表3所示，由表3可见，三个采样时刻的ASE均较小；利用步骤（4）判断得出预报警时刻k_b=278，温度的实际观测值在k_f=325时超出控制限，所以提前δk=k_f-k_b=47步预报，为维修人员采取有效措施排除故障提供了时间。

表3：平均归一化平方预测误差

Claims (7)

1.一种基于模糊自适应预测的动态过程故障预报方法，其特征在于实现动态过程故障预报的具体过程包括以下步骤：

（1）、先采集故障参数的实时测量数据，使用模糊自适应无迹卡尔曼预测方法计算故障参数后再确定过程状态；

（2）、根据步骤（1）计算的故障参数的标准差判断是否进行未来多步预测；若进行未来多步预测，则进行步骤（3）；否则，返回执行步骤（1）；

（3）、使用模糊自适应无迹卡尔曼预测方法计算未来多步预测值，利用未来测量输出变量的西格玛点计算故障预测概率；

（4）、根据故障预测概率判断动态过程将来是否超出正常运行区域；若超出正常运行区域，则进行故障预报警后返回执行步骤（1）；否则直接返回执行步骤（1）。

2.根据权利要求1所述的基于模糊自适应预测的动态过程故障预报方法，其特征在于涉及的模糊自适应无迹卡尔曼预测方法采用非线性动态系统模型，非线性动态系统模型为：

x(k+1)＝f(x(k),u(k))+d(k)+w(k)  （1）

y(k)＝H(k)x(k)+v(k)          （2）

其中k表示采样时刻，k为自然数，1≤k≤1000，x表示n维的状态向量，n为自然数，1≤n≤10，u表示l维的测量输入向量，l为自然数，1≤l≤10，d表示n维的不可测故障项，f(·)表示n维非线性函数向量，对x连续可导，y表示m维的测量输出向量，m为自然数，1≤m≤10，H表示m×n维的矩阵，w表示n维的过程噪声，服从均值为零协方差矩阵为Q的正态分布，Q表示n×n维的协方差矩阵，v表示m维的测量噪声，服从均值为零协方差矩阵为R的正态分布，R表示m×m维的协方差矩阵，w与v不相关；u(k)∈R^l×1，x(k)∈R^n×1，y(k)∈R^m×1分别表示测量输入、状态和测量输出；d(k)∈R^n×1表示不可测故障项；w(k)∈R^n×1，v(k)∈R^m×1分别表示过程噪声和测量噪声，二者不相关且分别服从均值为零、方差为Q(k)和R(k)的正态分布；当过程正常运行时，故障项d(k)＝0；当过程出现故障时，故障项d(k)≠0；故障预报针对缓慢变化型故障采用的故障模型为：

d(k)＝Γ(k,k₀)θ       （3）

其中θ表示2n维的故障参数向量，k₀表示故障起始时间，Γ(k,k₀)为n×2n维的时间函数矩阵；以2阶系统为例，故障模型为：

$d\left(k\right)=\mathrm{\&Gamma;}(k,k_0)\mathrm{\&theta;}=\left[\begin{array}{c}{\mathrm{\&theta;}}_{\mathrm{1,0}}+{\mathrm{\&theta;}}_{1,s}{T}_s(k-k_0)\\ {\mathrm{\&theta;}}_{\mathrm{2,0}}+{\mathrm{\&theta;}}_{2,s}{T}_s(k-k_0)\end{array}\right]---\left(4\right)$

其中

$\mathrm{\&Gamma;}(k,k_0)=\left[\begin{array}{cccc}1& 0& T_s(k-k_0)& 0\\ 0& 1& 0& T_s(k-k_0)\end{array}\right]$

θ＝[θ_1,0θ_2,0θ_1,sθ_2,s]^T

其中T_s为采样周期，0≤T_s≤1分钟，θ_h,0表示故障参数向量θ的第h个元素，θ_h,s表示故障参数向量θ的第2+h个元素，h＝1,2；

故障参数θ的动态模型为

$\left\{\begin{array}{c}\mathrm{\&theta;}(k+1)=\mathrm{\&theta;}\left(k\right)+\mathrm{\&epsiv;}\left(k\right)\\ y(k+1)=f(x\left(k\right),u\left(k\right))+\mathrm{\&Gamma;}(k,k_0)\mathrm{\&theta;}(k+1)+w^{\&prime;}\left(k\right)\end{array}\right.---\left(5\right)$

其中ε表示2n维的过程噪声，服从均值为零协方差矩阵为Q_θ的正态分布，Q_θ表示2n×2n维的协方差矩阵，w′(k)＝w(k)-Γ(k,k₀)ε(k)，w′表示m维的测量噪声，服从均值为零协方差矩阵为R_θ的正态分布，R_θ表示m×m维的协方差矩阵。

3.根据权利要求1所述的基于模糊自适应预测的动态过程故障预报方法，其特征在于涉及的模糊自适应无迹卡尔曼预测方法的工艺步骤包括单步预测、模糊自适应更新和未来多步预测三个步骤，其具体过程如下：

（A1）、单步预测：对公式（1）（2）所示的非线性动态系统模型，构造增广状态向量x_a(k-1)＝[x^T(k-1)  w^T(k-1)  v^T(k-1)]^T，设x_a(k-1)的第i个西格玛点为X_i(k-1)，x_a(k-1)的第i个西格玛点的权值为ω_i，X_i(k-1)为n_a维向量，0≤i≤2n_a，n_a＝2n+m，a表示增广，X_i(k-1)和ω_i的计算公式为：

$X_i(k-1)=\left\{\begin{array}{cc}{\tilde{x}}_a(k-1),& i=0\\ {\tilde{x}}_a(k-1)+A_i(k-1),& i=1,...,{2n}_a\\ {\tilde{x}}_a(k-1)-A_{i-n_a}(k-1),& i=n_a+1,...,{2n}_a\end{array}\right.{\mathrm{\&omega;}}_i=\left\{\begin{array}{cc}\frac{\mathrm{\&eta;}}{n_a+\mathrm{\&eta;}},& i=0\\ \frac{1}{2(n_a+\mathrm{\&eta;})},& i=1,...,{2n}_a\end{array}\right.---\left(6\right)$

其中 ${\tilde{x}}_a(k-1)={\left[\begin{array}{ccc}{\hat{x}}^T(k-1)& 0_{1\&times;n}& 0_{1\&times;m}\end{array}\right]}^T,$ 表示第k-1个采样时刻状态x的估计值，0_1×n表示1行n列的零矩阵，0_1×m表示1行m列的零矩阵，向量A_i表示矩阵A的第i列，矩阵A为n_a×n_a维，满足A(k-1)A^T(k-1)＝(n_a+η)P_a(k-1)，矩阵P_a(k-1)为n_a×n_a维，表示第k-1个采样时刻增广状态的误差协方差矩阵，参数η根据状态x的分布确定，对于高斯分布η＝3-n_a；记预测状态的第i个西格玛点为

状态预测值为其误差协方差矩阵为P(k|k-1)，其中

为n维向量，

表示在第k-1个采样时刻预测的第k个采样时刻的x的估计值，P(k|k-1)为n×n维矩阵，表示在第k-1个采样时刻预测的第k个采样时刻的误差协方差矩阵，则

和P(k|k-1)的计算公式如下：

${\tilde{X}}_i\left(k\right)=f({X_i}^x(k-1),u(k-1))+{X_i}^w(k-1)+\mathrm{\&Gamma;}(k-1,k_0)\widehat{\mathrm{\&theta;}}(k-1)---\left(7\right)$

$\hat{x}\left(k\right|k-1)=\underset{i=0}{\overset{{2n}_a}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_i{\tilde{X}}_i\left(k\right)---\left(8\right)$

$P\left(k\right|k-1)=\underset{i=0}{\overset{{2n}_a}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_i[{\tilde{X}}_i\left(k\right)-\hat{x}\left(k\right|k-1)]{[{\tilde{X}}_i\left(k\right)-\hat{x}\left(k\right|k-1)]}^T---\left(9\right)$

其中

为故障参数θ(k-1)的估计值，

和

分别表示X_i(k-1)中与状态和过程噪声对应的部分，

由Xi(k-1)的第1至n个元素组成，由X_i(k-1)的第n+1至2n个元素组成；记预测输出的第i个西格玛点为

输出预测值为

其中

为m维向量，

表示在第k-1个采样时刻预测的第k个采样时刻的y的估计值，则

的计算公式为：

${\tilde{Y}}_i\left(k\right|k-1)=H\left(k\right){\tilde{X}}_i\left(k\right)+{X_i}^v(k-1)---\left(10\right)$

$\hat{y}\left(k\right|k-1)=\underset{i=0}{\overset{{2n}_a}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_i{\tilde{Y}}_i\left(k\right|k-1)---\left(11\right)$

其中X_i ^v(k-1)表示X_i(k-1)中与测量噪声对应的部分，X_i ^v(k-1)由X_i(k-1)的第2n+1至2n+m个元素组成；

（A2）、模糊自适应更新：利用输出预测误差

及其协方差矩阵的估计值

计算自适应调节参数λ(k)，

$P_{\tilde{y}}\left(k\right)=\left\{\begin{array}{cc}V\left(k\right),& k=1\\ \frac{\mathrm{\&rho;V}(k-1)+V\left(k\right)}{1+\mathrm{\&rho;}},& k\&GreaterEqual;2\end{array}\right.---\left(12\right)$

$\mathrm{\&lambda;}\left(k\right)=\max \{\frac{\mathrm{tr}[P_{\tilde{y}}\left(k\right)-\mathrm{\&zeta;R}\left(k\right)]}{\mathrm{tr}\left[H\left(k\right)P\left(k\right|k-1)H^T\left(k\right)\right]},1\}---\left(13\right)$

其中

ρ为遗忘因子，0≤ρ≤1，ζ为弱化因子，1≤ζ≤10，tr[·]为求迹算子；由于过程故障运行状态的模糊性，利用模糊逻辑系统计算自适应调节参数λ^*(k)，令模糊集合“Normal”和“Failure”分别表示正常运行状态和故障运行状态，对应的隶属度函数分别为μ_N(y(k))和μ_F(y(k))，且μ_N(y(k))+μ_F(y(k))＝1；模糊规则为：

R¹:IFy(k)is“Normal”,THEN

R²:IFy(k)is“Failure”,THEN

设测量输出向量y(k)＝[y₁(k),…y_z(k),…y_m(k)]^T服从均值为μ协方差为Σ的正态分布，其中y_z(k)为y(k)的第z个元素，1≤z≤m，μ为m维向量，Σ为m×m维矩阵，则

$\begin{array}{c}{\mathrm{\&mu;}}_F\left(y\left(k\right)\right)=P(-|y_1\left(k\right)|\&le;{\mathrm{\&psi;}}_1<|y_1\left(k\right)|,\&CenterDot;\&CenterDot;\&CenterDot;-|y_z\left(k\right)|\&le;{\mathrm{\&psi;}}_z<|y_z\left(k\right)|,\&CenterDot;\&CenterDot;\&CenterDot;,-|y_m\left(k\right)|\&le;{\mathrm{\&psi;}}_m<|y_m\left(k\right)\left|\right)\\ =\frac{1}{\sqrt{{\left(2\mathrm{\&pi;}\right)}^m\left|\mathrm{\&Sigma;}\right|}}{\&Integral;}_{-\left|y_1\left(k\right)\right|}^{\left|y_1\left(k\right)\right|}\&CenterDot;\&CenterDot;\&CenterDot;{\&Integral;}_{-\left|y_z\left(k\right)\right|}^{\left|y_z\left(k\right)\right|}\&CenterDot;\&CenterDot;\&CenterDot;{\&Integral;}_{-\left|y_m\left(k\right)\right|}^{\left|y_m\left(k\right)\right|}{e}^{-\frac{1}{2}{(\mathrm{\&psi;}-\mathrm{\&mu;})}^T{\mathrm{\&Sigma;}}^{-1}(\mathrm{\&psi;}-\mathrm{\&mu;})}\mathrm{d\&psi;}\end{array}---\left(14\right)$

其中向量ψ＝[ψ₁,…ψ_z,…,ψ_m]^T，ψ_z为ψ的第z个元素，P(·)表示概率函数，利用Monte-Carlo方法或分区自适应方法等数值积分算法计算式（14）；先计算自适应调节参数λ^*(k)

${\mathrm{\&lambda;}}^{*}\left(k\right)={\mathrm{\&mu;}}_N\left(y\left(k\right)\right){\mathrm{\&lambda;}}_1^{*}\left(k\right)+{\mathrm{\&mu;}}_F\left(y\left(k\right)\right){\mathrm{\&lambda;}}_2^{*}\left(k\right)---\left(15\right)$

再利用λ^*(k)调节P(k|k-1)

$\stackrel{\&OverBar;}{P}\left(k\right|k-1)={\mathrm{\&lambda;}}^{*}\left(k\right)P\left(k\right|k-1)---\left(16\right)$

其中矩阵表示矩阵P(k|k-1)的模糊自适应修正值，则卡尔曼增益K(k)、状态修正值及其误差协方差矩阵P(k)的计算公式为：

$K\left(k\right)=\stackrel{\&OverBar;}{P}\left(k\right|k-1)H^T\left(k\right){[H\left(k\right)\stackrel{\&OverBar;}{P}\left(k\right|k-1)H^T\left(k\right)+R\left(k\right)]}^{-1}---\left(17\right)$

$\hat{x}\left(k\right)=\hat{x}\left(k\right|k-1)+K\left(k\right)(y\left(k\right)-\hat{y}\left(k\right|k-1))---\left(18\right)$

$P\left(k\right)=\stackrel{\&OverBar;}{P}\left(k\right|k-1)-K\left(k\right)[H\left(k\right)\stackrel{\&OverBar;}{P}\left(k\right|k-1)H^T\left(k\right)+R\left(k\right)]K^T\left(k\right)---\left(19\right);$

（A3）、未来多步预测；在第k个采样时刻，重复迭代执行步骤（A1）得到未来p步的输出预测值

表示在第k个采样时刻预测的第k+j个采样时刻的y的预测值，j＝1,2,…,p，p为未来预测步数，1≤p≤100，其中

P(k+j)＝P(k+j|k)（20）

$\hat{x}(k+j)=\hat{x}(k+j|k)---\left(21\right)$

${\tilde{Y}}_i(k+j|k)=H(k+j){\tilde{X}}_i(k+j)+{X_i}^v(k+j-1)---\left(22\right)$

P(k+j|k)表示在第k个采样时刻预测的第k+j个采样时刻的误差协方差矩阵，

表示在第k个采样时刻预测的第k+j个采样时刻的x的估计值，

表示在第k个采样时刻预测的第k+j个采样时刻的西格玛点。

4.根据权利要求1所述的基于模糊自适应预测的动态过程故障预报方法，其特征在于涉及的步骤（1）确定过程状态时，在每个采样时刻，先计算故障参数的估计值

然后确定状态x，其具体过程为：

（1.1）、在第k个采样时刻，根据故障参数的动态模型公式（5），将故障参数作为状态，利用模糊自适应无迹卡尔曼预测方法的步骤（A1）和步骤（A2）计算故障参数的估计值

或利用扩展卡尔曼滤波器估计；

（1.2）、根据过程动态模型公式（1）和（2），利用模糊自适应无迹卡尔曼预测方法的步骤（A1）和步骤（A2）计算过程状态的估计值

5.根据权利要求1所述的基于模糊自适应预测的动态过程故障预报方法，其特征在于涉及的步骤（2）中根据故障参数的标准差判断是否进行未来多步预测，其判断准则为：在第k个采样时刻，若 $\mathrm{std}\left({\widehat{\mathrm{\&theta;}}}_{\mathrm{1,0}}\left(k\right)\right)<{\mathrm{\&sigma;}}_1$ 且…且 $\mathrm{std}\left({\widehat{\mathrm{\&theta;}}}_{\mathrm{1,0}}\left(k\right)\right)<{\mathrm{\&sigma;}}_h$ 且…且 $\mathrm{std}\left({\widehat{\mathrm{\&theta;}}}_{\mathrm{1,0}}\left(k\right)\right)<{\mathrm{\&sigma;}}_n,$ 其中std(·)表示标准差函数，σ_h表示标准差的阈值，根据过程噪声w的标准差确定，h＝1,2,…,n，则进行未来预测的采样时刻k_p＝k，执行步骤（3）计算未来p步预测值；否则，返回执行步骤（1），继续跟踪过程的变化趋势；其中阈值σ_h根据过程噪声w的标准差确定。

6.根据权利要求1所述的基于模糊自适应预测的动态过程故障预报方法，其特征在于涉及的步骤（3）中未来多步预测值和故障预测概率的计算方法为：

（3.1）、在第k_p个采样时刻，利用模糊自适应无迹卡尔曼预测方法的步骤（A3）计算未来p步的输出预测值

及其西格玛其中j＝1,2,…,p，i＝0,1,2,…,2n_a；

（3.2）、对未来输出预测值的西格玛进行分析，若

则则第i个故障指标函数F_i(k_p+j|k_p)＝1；否则F_i(k_p+j|k_p)＝0；其中

为第i个西格玛的第z个元素；B_z为第z个输出变量的故障区间，由输出变量的控制限决定；

i＝0,1,2,…,2n_a,z＝1,2,…,m；按照公式（24）计算故障预测概率：

$P_f(k_p+j|k_p)=\underset{i=0}{\overset{{2n}_a}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_i{F}_i(k_p+j|k_p)---\left(24\right).$

7.根据权利要求1所述的基于模糊自适应预测的动态过程故障预报方法，其特征在于涉及的步骤（4）中故障预报警判断准则为：在第k_p个采样时刻，根据P_f(k_p+j|k_p)判断过程将来是否超出正常运行区域发生故障，若P_f(k_p+q|k_p)＝1且P_f(k_p+q+g|k)＝1，则预报警采样时刻k_b＝k_p，发出故障预报警，然后返回执行步骤（1）；否则，直接返回执行步骤（1），继续跟踪过程变化趋势；其中q＝1,2,…,p-1,g＝1,2,…,p-q。

Images