CN105676833A - 发电过程控制系统故障检测方法 - Google Patents

Abstract

本发明公开了一种发电过程控制系统故障检测方法，所述方法包括：用主成分分析PCA法对经过降噪、标准化处理后的训练矩阵X进行矩阵分解，采用得分矩阵T作为基矩阵W的初始值W₀；用带非负约束的交替最小二乘法迭代求解所述训练矩阵X的基矩阵W和权重系数矩阵H；构造基于非负矩阵分解的监控统计量和SPE_n，利用核密度估计法分别计算监控统计量和SPE_n的概率密度函数PDF，设置显著水平并分别求取统计量和SPE_n的控制限；利用所述权重系数矩阵H和经过数据处理之后的测试矩阵X_test计算得出所述测试矩阵X_test基矩阵的近似值分别计算所述测试矩阵X_test的监控统计量和SPE_n，并与相应的所述控制限比较，超过所述控制限则表明有故障发生。本发明能够针对发电过程海量运行数据进行状态监测，进而实现发电过程控制系统的故障诊断。

Description

发电过程控制系统故障检测方法

技术领域

本发明涉及发电技术领域，尤其涉及一种发电过程控制系统故障检测方法。

背景技术

发电过程是一类典型的复杂工业过程，对其控制系统进行故障检测与诊断已经成为控制理论研究的一个重要分支。对于复杂工业过程来说想要建立其精确的数学模型是非常困难的，但是发电过程分散控制系统(Distributedcontrolsystem，DCS)与监控信息系统(Supervisoryinformationsystem，SIS)已实现了生产过程的数字化，使运行数据得以长期保存。因此，利用机组正常运行的历史数据建立监视系统运行状态的模型是一个既方便又有效的方法。

如何从大量、高维的运行数据中提取有效的故障特征信息将是至关重要的问题。矩阵分解技术广泛应用于各类学科的研究中，通过矩阵分解可以将原始数据从高维空间投影到低维空间，从而发现原始数据的内在结构特征。常见的传统的矩阵分解方法有：主成分分析(Principalcomponentanalysis，PCA)，快速独立主成分分析(FastIndependentcomponentanalysis，FastICA)，Fisher判据分析(Fisherdiscriminantanalysis，FDA)等。在这些方法的计算过程中没有对数据进行非负约束使得其计算结果当中可能包含负数，这在数值计算的角度来看是没有问题的。但是，对于绝大多数工业过程来说，其运行数据都是非负的，这就在一定程度上造成了计算结果的不可解释性。同时，传统的PCA方法假设过程变量是服从高斯分布的，然而在复杂工业过程中这个假设很难成立的；FastICA算法对初始值的选择要求较高，如果初始值选择的不合适有可能会造成算法不收敛；FDA方法在实现过程中对数据均值信息的依赖程度很高，当处理均值变化不大的数据集时效果不理想。

为此，需要一种实现简便、分解形式和分解结果具有可解释性、以及占用存储空间少的发电过程控制系统故障检测方法。

发明内容

为解决现有存在的技术问题，本发明实施例提供一种发电过程控制系统故障检测方法。

为达到上述目的，本发明实施例的技术方案是这样实现的：

一种发电过程控制系统故障检测方法，所述方法包括：

用主成分分析PCA法对经过降噪、标准化处理后的训练矩阵X进行矩阵分解，采用得分矩阵T作为基矩阵W的初始值W₀；

用带非负约束的交替最小二乘法迭代求解所述训练矩阵X的基矩阵W和权重系数矩阵H；

构造基于非负矩阵分解的监控统计量和SPE_n，利用核密度估计法分别计算监控统计量和SPE_n的概率密度函数PDF，设置显著水平并分别求取统计量和SPE_n的控制限；

利用所述权重系数矩阵H和经过数据处理之后的测试矩阵X_test计算得出所述测试矩阵X_test基矩阵的近似值分别计算所述测试矩阵X_test的监控统计量和SPE_n，并与相应的所述控制限比较，超过所述控制限则表明有故障发生。

其中，所述用主成分分析PCA法对经过降噪、标准化处理后的训练矩阵X进行矩阵分解，采用得分矩阵T作为基矩阵W的初始值W₀，包括：

采集一个系统正常运行时的样本集合矩阵作为所述训练矩阵X，X∈R^n×m，n为数据样本个数，m为变量个数；

对所述训练矩阵X进行数据预处理：用带有遗忘因子的递推方法对数据样本进行降噪；对降噪后的所述训练矩阵X做标准化处理；

采用PCA方法对所述训练矩阵X进行矩阵分解；

采用基于故障信噪比确定所述得分矩阵T的主元个数a，并求得所述得分矩阵T，即基矩阵W的初始值W₀。

其中，所述用带非负约束的交替最小二乘法迭代求解所述训练矩阵X的基矩阵W和权重系数矩阵H，包括：将主元分解得到的基矩阵初始值W₀标准化；

在传统的交替最小二乘法中引入两个稀疏因子，通过增加所述两个稀疏因子的提高所述基矩阵W和权重系数矩阵H的稀疏程度；运用带约束的交替最小二乘法求解所述基矩阵W和权重系数矩阵H。

其中，所述构造基于非负矩阵分解的监控统计量和SPE_n，利用核密度估计法分别计算监控统计量和SPE_n的概率密度函数，包括：对所述基矩阵W进行重构，定义重构后的基矩阵为：其中，H为权重系数矩阵，T为得分矩阵，λ_W为稀疏因子；基于NMF的监控模型将所述训练矩阵X描述为E为残差矩阵；定义基于NMF的监控统计量和SPE_n：其中，表示一个样本向量的重构值，II表示a×a的单位矩阵，a为所述得分矩阵T的主元个数；

采用所述核密度估计法估计所述监控统计量和SPE_n的PDF，采用高斯核函数为核函数，带宽由广义交叉熵算法求取。

采用所述核密度估计法估计所述监控统计量和SPE_n的PDF，包括：

定义训练矩阵X的密度函数f(x)的估计形式为：其中，n是样本个数，h是带宽，K(·)是核函数，满足：u表示变量，x表示需要求概率密度函数的变量的数值，x_i表示数列元素；

使下列Csiszár测度达到最小值，即：其中，p(x)为给定数据的先验分布概率，若先验概率未知，则p(x)＝1，g(x)是的另一种表现形式，λ＝[λ₁,…,λ_n]^T为拉格朗日乘子；

根据求解得到广义交叉熵问题中带宽的解h^*，其中，C为n×n方阵，其元素为熵值，

C和均为带宽h和数据x的函数；

通过将h^*带入式计算h_opt；

根据h_opt以及式计算得到密度函数f(x)，密度函数f(x)即为要估算的PDF。

本发明实施例提供一种发电过程控制系统故障检测方法，能够针对发电过程海量运行数据进行状态监测，进而实现发电过程控制系统的故障诊断，为控制系统优化设计与调试维护提供了先进的检测工具。

附图说明

在附图(其不一定是按比例绘制的)中，相似的附图标记可在不同的视图中描述相似的部件。具有不同字母后缀的相似附图标记可表示相似部件的不同示例。附图以示例而非限制的方式大体示出了本文中所讨论的各个实施例。

图1为本发明实施例发电过程控制系统故障检测方法的流程图。

具体实施方式

本发明提供一种基于非负矩阵分解的发电过程控制系统故障检测方法，采用OPC通讯的方式与SIS连接，从SIS数据库中获取发电过程实时或历史运行数据或者在SIS站中进行二次开发也可以实现该方法的所有功能，其流程图如图1所示，包括以下步骤：

步骤101、用PCA对经过降噪、标准化处理后的训练矩阵X进行矩阵分解，采用得分矩阵T作为基矩阵W的初始值W₀；

步骤102、用带非负约束的交替最小二乘法(ALS，AlternatingLeast-squares)迭代求解X的基矩阵W和权重系数矩阵H；

步骤103、构造基于非负矩阵分解的监控统计量和SPE_n，利用核密度估计法(KDE，KernelDensityEstimation)分别计算监控统计量和SPE_n的概率密度函数(PDF，ProbabilityDensityFunction)，设置显著水平并分别求取统计量和SPE_n的控制限；

步骤104、利用步骤102中分解得到权重系数矩阵H和经过数据处理之后的测试矩阵X_test计算得出X_test基矩阵的近似值分别计算X_test的监控统计量和SPE_n并与相应的控制限进行比较，如果超过控制限，则表明有故障发生，反之表明系统运行正常。

其中，步骤101具体可以包括：

步骤11：采集一个系统正常运行时的样本集合矩阵作为训练矩阵X(X∈R^n×m)，n为数据样本个数，m为变量个数；

步骤12：对训练矩阵X进行数据预处理；

1)用带有遗忘因子的递推方法对数据样本进行降噪，其过程可用下式(1)来描述：

x(k)＝λx(k)+(1-λ)x(k-1)(1)

其中，x(k)为第k时刻变量的测量值，λ为遗忘因子，0≤λ≤1，一般取0.4～0.5。

2)对降噪后的训练矩阵X做标准化处理，计算公式如下式(2)：

$x\left(j\right)=\frac{x\left(j\right)-\stackrel{\‾}{x}\left(j\right)}{s\left(j\right)}---\left(2\right)$

其中，x(j)为训练集的列向量，为变量均值，s(j)为变量的标准差，s(j)的计算如下式(3)所示：

$s\left(j\right)=\sqrt{\frac{\underset{i=1}{\overset{n}{\Σ}}{(x_{ij}-{\stackrel{\‾}{x}}_j)}^2}{n-1}}---\left(3\right)$

步骤13：采用PCA方法对训练矩阵X进行矩阵分解；

定义标准化后的样本集合矩阵X协方差矩阵为下式(4)：

$S=\mathrm{cov}\left(X\right)\≈\frac{1}{n-1}{X}^{T}X---\left(4\right)$

对协方差矩阵S进行奇异值分解，S＝U∑V^T。

采用PCA对训练矩阵X进行如下分解：

$X=\hat{X}+E={\mathrm{TP}}^T+E---\left(5\right)$

T＝XP(6)

其中，P∈R^m×a为负载矩阵，由V的前a个奇异向量构成，T∈R^n×a为得分矩阵，T的各列被称为主元变量，a表示主元个数，也是得分矩阵的列数，E为残差矩阵，U是进行奇异值分解之后得到的左奇异矩阵；

步骤14：采用基于故障信噪比确定主元个数a；

将某一被测量发生的故障记为fξ_i，则被测量的输出可以描述为下式(7)：

x＝x^*+fξ_i(7)

其中，x^*是被测量正常时的观测值。f是一个标量值，表示故障的幅值。ξ_i是故障方向向量，描述故障的程度。当被测量发生故障时SPE统计量可以写成下式(8)：

$SPE={||x^T-x^T{P}_a{P}_a^T||}^2={||x^T{\tilde{P}}_a{\tilde{P}}_a^T||}^2---\left(8\right)$

其中，P_a表示主元个数为a时的负载矩阵，将式(7)带入式(8)，可以得到下式(9)：

$\begin{array}{c}SPE={||{(x^{*}+{\mathrm{f\ξ}}_i)}^T-{(x^{*}+{\mathrm{f\ξ}}_i)}^T{P}_a{P}_a^T||}^2\\ ={||{(x^{*}+{\mathrm{f\ξ}}_i)}^T{\tilde{P}}_a{\tilde{P}}_a^T||}^2={||{\mathrm{\ξ}}_i^T{\tilde{P}}_a{\tilde{P}}_a^T||}^2\end{array}---\left(8\right)$

由于x^*标准化后其均值为0、方差为1，因此，令x^*＝0，f＝1，定义SPE统计量的故障信噪比为下式(10)：

${\mathrm{SNR}}_Q=\frac{{||{\mathrm{\ξ}}_i^T{\tilde{P}}_a{\tilde{P}}_a^T||}^2}{Q_{\mathrm{\α}}}=\frac{{||{\widetilde{\mathrm{\ξ}}}_i||}^2}{Q_{\mathrm{\α}}}---\left(10\right)$

其中，是故障方向ξ_i在残差空间投影的平方范数，Q_α是Q统计量的控制限。

同理，将式(7)带入T²＝x^TPΛ^-1P^Tx得到T²统计量的描述为下式(11)：

T²＝(x^*+fξ_i)^TP_aλ^-1P_a ^T(x^*+fξ_i)(11)

同样，定义T²统计量的故障信噪比为下式(12)：

${\mathrm{SNR}}_{T^2}=\frac{{||{(x^{*}+{\mathrm{f\ξ}}_i)}^T{P}_a{\mathrm{\λ}}^{-1}{P_a}^T(x^{*}+{\mathrm{f\ξ}}_i)||}^2}{T_{a,m,\mathrm{\α}}^2}=\frac{{||{\mathrm{\ξ}}_i^T{P}_a{\mathrm{\λ}}^{-1}{P_a}^T{\mathrm{\ξ}}_i||}^2}{T_{\mathrm{\α},m,\mathrm{\α}}^2}---\left(12\right)$

其中，分子是故障方向的T²统计量信息，分母是T²统计量的控制限。

从故障检测的几何意义上来讲，统计量与其相应的控制限的比值反映了故障检测的灵敏度，因此，使得故障信噪比取得最大值的主元个数，即为最优主元个数。由T＝XP便可求得得分矩阵，即基矩阵W的初始值W₀。

其中，步骤102具体可以包括：

步骤21：将主元分解得到的基矩阵初始值W₀标准化；

步骤22：在传统的交替最小二乘法中引入两个稀疏因子λ_H和λ_W(λ_H＞0，λ_W＞0)，通过增加λ_H和λ_W的值可以提高基矩阵W和权重系数矩阵H的稀疏程度。增加约束之后的最小二乘问题可以描述为下式(13)：

$\begin{array}{ccc}\underset{h_j}{\min }{||x_j-{\mathrm{Wh}}_j||}_2^2+{\mathrm{\λ}}_H{||h_j||}_2^2& s.t.& {\mathrm{\λ}}_H\≥0,h_j\≥0\end{array}---\left(13\right)$

其中，x_j和h_j分别代表X和H的列向量。

步骤23：运用带约束的交替最小二乘法求解基矩阵W和权重系数矩阵H。

1)设定稀疏因子λ_H和λ_W的值，以及算法最大迭代次数maxiter；

2)固定W，将W₀带入(W^TW+λ_HI)H＝W^TX中求解权重系数矩阵H，I表示a×a的单位矩阵，a为所述所述得分矩阵T的主元个数；

3)检查矩阵H中所有元素的大小，将小于0的元素设为0；

4)固定H，将3)中得到的H带入(HH^T+λ_WI)W^T＝HX^T中求解基矩阵W；

5)检查矩阵W中所有元素的大小，将小于0的元素设为0；

6)判断算法迭代次数是否达到最大次数，是跳出循环，不是继续执行下一步；

7)计算W每列的L₂范数，将W的列向量按其L₂范数的大小降序排列；

8)将矩阵W标准化；

9)重复执行步骤2)—8)直到算法跳出循环，输出W和H的最终值。

其中，步骤103具体可以包括：

步骤31：对基矩阵W进行重构，定义重构后的基矩阵为下式(14)：

$\hat{W}=X{\left({\left({\mathrm{HH}}^T+{\mathrm{\λ}}_{W}I\right)}^{-1}H\right)}^T---\left(14\right)$

步骤32：基于NMF的监控模型描述训练矩阵X为下式(15)：

$X=\hat{W}H+E---\left(15\right)$

参考基于主元分析故障诊断方法中对监控统计量T²和SPE的定义，定义基于NMF的监控统计量和SPE_n为下式(16)、(17)所示：

$T_n^2\left(i\right)=\hat{W}\left(i\right){\mathrm{\Λ}}^{-1}\hat{W}{\left(i\right)}^T---\left(16\right)$

${\mathrm{SPE}}_n\left(i\right)=\left(x\right(i)-\hat{x}(i\left)\right){\left(x\right(i)-\hat{x}(i\left)\right)}^T---\left(17\right)$

其中，表示一个样本向量的重构值，计算如下式(18)：

$\hat{x}\left(i\right)=\hat{W}\left(i\right)H=x\left(i\right){\left({\left({\mathrm{HH}}^T+{\mathrm{\λ}}_{W}I\right)}^{-1}H\right)}^{T}H---\left(18\right)$

步骤33：采用KDE方法估计监控统计量和SPE_n的概率密度函数(PDF)。在现实中，训练矩阵X的密度函数f(x)是不能准确得到的，定义其估计形式如下式(19)：

$\hat{f}\left(x\right)=\frac{1}{nh}\underset{i=1}{\overset{n}{\Σ}}K\left(\frac{x_i-x}{h}\right)---\left(19\right)$

其中，n是样本个数，h是带宽，K(·)是核函数，且满足下式(20)：

$\begin{array}{cc}{\∫}_{-\mathrm{\∞}}^{\mathrm{\∞}}K\left(u\right)du=1& K\left(u\right)\≥0\end{array}---\left(20\right)$

在KDE方法中核函数和带宽是需要确定的参数，本申请采用高斯核函数为核函数，带宽由广义交叉熵(GCE)算法求取。

GCE算法的求解目标是使下列Csiszár测度达到最小值，如下式(21)：

$D(g\→p)=-\frac{1}{2}+\frac{1}{2}\∫\frac{g^2\left(x\right)}{p\left(x\right)}dx---\left(21\right)$

其中，Csiszár测度可以理解成表示一种距离，D(g→p)即为Csiszár测度，p(x)为给定数据的先验分布概率，若先验概率未知，则p(x)＝1。g(x)是的另一种表现形式，如下式(22)：

$g\left(x\right)=p\left(x\right)\underset{j=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_{j}K(x;x_j,h)---\left(22\right)$

其中，λ＝[λ₁,…,λ_n]^T为拉格朗日乘子。

求解GCE问题即求解拉格朗日乘子λ和带宽h，通过凸二次规划问题(CQPP)给出求解式(21)最小化的转化形式如下式(23)：

$(h^{*},{\mathrm{\λ}}^{*})=\left\{\right(h,\mathrm{\λ}):c^T\mathrm{\λ}(h)=1,\mathrm{\λ}(h)=\underset{\mathrm{\λ}\≥0}{\arg \min }(\frac{1}{2}{\mathrm{\λ}}^{T}C\left(h\right)\mathrm{\λ}-{\mathrm{\λ}}^T\hat{k}\left(h\right)\left)\right\}---\left(23\right)$

其中，C为n×n方阵，其元素为熵值，分别计算如下式(24)：

${\hat{k}}_i\left(h\right)=\frac{1}{n-1}\underset{j\≠i}{\overset{n}{\Σ}}K(x_j;h,x_i)---\left(25\right)$

由上式可知，C和均为带宽h和数据矩阵x的函数，u表示变量，x表示需要求概率密度函数的变量的数值，x_i表示数列元素，则通过求解式(23)可以得到GCE问题中带宽的解h^*，但这并不是KDE问题中式(19)的最优带宽值，要求得式(19)中的带宽最优值还要将h^*带入下式(26)计算h_opt：

$h_{opt}={\left(2n\sqrt{\mathrm{\π}}{||{\left(g\right(\·;h^{*})/p)}^{\′\′}||}^2\right)}^{-2/5}---\left(26\right)$

此时，根据h_opt以及式(19)计算得到密度函数f(x)，密度函数f(x)即为要估算的PDF，如此，将上文计算得到的监控统计量和SPE_n的值带入KDE方法中便可得到其分布情况。设置显著性水平为0.99(实际应用中，显著水平的设置取值还可以是0.95或其他，可以通过查表得到)，分别计算监控统计量和SPE_n的控制限和SPE_nα。

其中，步骤104具体可以包括：

根据上文求解得到的权重系数矩阵H和测试矩阵经过数据处理之后的X_test利用式(14)计算得出X_test基矩阵的近似值分别计算X_test的监控统计量和SPE_n并与相应的控制限进行比较，如果超过控制限则表明有故障发生，反之表明系统运行正常。

以上所述，仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。

Claims (5)

1.一种发电过程控制系统故障检测方法，其特征在于，所述方法包括：

用主成分分析PCA法对经过降噪、标准化处理后的训练矩阵X进行矩阵分解，采用得分矩阵T作为基矩阵W的初始值W₀；

用带非负约束的交替最小二乘法迭代求解所述训练矩阵X的基矩阵W和权重系数矩阵H；

构造基于非负矩阵分解的监控统计量和SPE_n，利用核密度估计法分别计算监控统计量和SPE_n的概率密度函数PDF，设置显著水平并分别求取统计量和SPE_n的控制限；

利用所述权重系数矩阵H和经过数据处理之后的测试矩阵X_test计算得出所述测试矩阵X_test基矩阵的近似值分别计算所述测试矩阵X_test的监控统计量和SPE_n，并与相应的所述控制限比较，超过所述控制限则表明有故障发生。

2.根据权利要求1所述的方法，其特征在于，所述用主成分分析PCA法对经过降噪、标准化处理后的训练矩阵X进行矩阵分解，采用得分矩阵T作为基矩阵W的初始值W₀，包括：

采集一个系统正常运行时的样本集合矩阵作为所述训练矩阵X，X∈R^n×m，n为数据样本个数，m为变量个数；

对所述训练矩阵X进行数据预处理：用带有遗忘因子的递推方法对数据样本进行降噪；对降噪后的所述训练矩阵X做标准化处理；

采用PCA方法对所述训练矩阵X进行矩阵分解；

采用基于故障信噪比确定所述得分矩阵T的主元个数a，并求得所述得分矩阵T，即基矩阵W的初始值W₀。

3.根据权利要求1所述的方法，其特征在于，所述用带非负约束的交替最小二乘法迭代求解所述训练矩阵X的基矩阵W和权重系数矩阵H，包括：

将主元分解得到的基矩阵初始值W₀标准化；

在传统的交替最小二乘法中引入两个稀疏因子，通过增加所述两个稀疏因子的提高所述基矩阵W和权重系数矩阵H的稀疏程度；

运用带约束的交替最小二乘法求解所述基矩阵W和权重系数矩阵H。

4.根据权利要求1所述的方法，其特征在于，所述构造基于非负矩阵分解的监控统计量和SPE_n，利用核密度估计法分别计算监控统计量和SPE_n的概率密度函数，包括：

对所述基矩阵W进行重构，定义重构后的基矩阵为：其中，H为权重系数矩阵，T为得分矩阵，λ_W为稀疏因子；

基于NMF的监控模型将所述训练矩阵X描述为E为残差矩阵；

定义基于NMF的监控统计量和SPE_n： 其中，表示一个样本向量的重构值，II表示a×a的单位矩阵，a为所述得分矩阵T的主元个数；

采用所述核密度估计法估计所述监控统计量和SPE_n的PDF，采用高斯核函数为核函数，带宽由广义交叉熵算法求取。

5.根据权利要求4所述的方法，其特征在于，所述采用所述核密度估计法估计所述监控统计量和SPE_n的PDF，包括：

定义训练矩阵X的密度函数f(x)的估计形式为：其中，n是样本个数，h是带宽，K(·)是核函数，满足：u表示变量，x表示需要求概率密度函数的变量的数值，x_i表示数列元素；

使下列Csiszár测度达到最小值，即：其中，p(x)为给定数据的先验分布概率，若先验概率未知，则p(x)＝1，g(x)是的另一种表现形式，λ＝[λ₁,…,λ_n]^T为拉格朗日乘子；

根据求解得到广义交叉熵问题中带宽的解h^*，其中，C为n×n方阵，其元素为熵值，

C和均为带宽h和数据x的函数；

通过将h^*带入式计算h_opt：；

根据h_opt：以及式计算得到密度函数f(x)，密度函数f(x)即为要估算的PDF。