CN104595170A - 一种自适应核高斯混合模型的空压机监控诊断系统及方法 - Google Patents

Abstract

本发明公开了一种自适应核高斯混合模型的空压机监控诊断系统及方法，涉及空压机控制技术领域。该系统包括现场设备层、设备控制层和管理监控层，现场设备层由PLC200、传感器、空压机、执行器和水泵构成，以PLC200作为从站，完成对现场设备层的控制；设备控制层包括上位机和PLC300，以PLC300作为主站，采用基于支持向量机的变结构自适应PID控制器对整个空压机系统进行控制，上位机实现对空压机系统的监控；设备控制层和管理监控层通过工业以太网进行通信，实现上位机的远程监控和数据传输。上位机上采用的故障诊断方法融合了高斯混合模型和核主成元分析，通过迭代优化方法寻找最优核函数参数，达到分开不同模态数据的目的。本发明具有更高的诊断精度和实用价值。

Description

一种自适应核高斯混合模型的空压机监控诊断系统及方法

技术领域

本发明涉及空压机控制技术领域，具体是一种自适应核高斯混合模型的空压机监控诊断系统及方法。

背景技术

空压机是煤矿重要的大型设备，也是高耗能设备。随着煤矿新工艺、动力设备及新技术的应用，不仅需要更多的风量，而且要求空压机能够随着负荷的变化而自适应地调节。而现有的老式空压机存在控制品质差、可靠性低以及难以实时监测系统工作状态等问题，造成了巨大的能源浪费，也给操作人员带来了巨大负担。传统的控制系统采用就地分散式人工操作。空压机工作时噪音很大，长时间工作在这种环境给工人的健康带来很大的危害，不够人性化。空压机每天24小时不间断运行需值班人员连续值守，控制水平的低下造成人力资源严重浪费，并且空压机需要连续运行，维护周期较长，空压机工作性能不稳定，易出现报警、停机等的误操作。

对煤矿空压机来说，目前人们还只停留在故障发生后的定性分析上，而对空压机的实时监测和故障诊断还比较少。绝大部分空压机监控系统只是对工艺变量进行上下限检查，以此判断空压机的工作状态。实际采集的数据含有大量噪声，这些噪声难以使用现有的数学模型描述，这类方法难以确定合适的阈值。阈值过大会导致难以准确监测轻微的故障，阈值太小会导致误报率过高，其次，由于空压机工作过程变量之间的很强的相关性，且变量之间的相关性随着工况的变化而变化，个别变量超出阈值并一定发生故障。再次，个别变量在设定的阈值范围内并一定表明系统没有发生故障，如传感器故障失效。近年来，作为基于知识的多变量统计控制(MSPC)得到了较为广泛的发展，而在MSPC中PCA方法是一种处理工业过程控制监测变量间的非自相关性的有效方法。但是PCA方法不能监测到输入数据的非线性结构。针对空压机传统的监测方法监测准确率低，误报率高等缺点，提出一种自适应核高斯混合模型的空压机监控诊断系统及方法。

发明内容

为了克服上述现有技术的缺点，本发明的一种一种自适应核高斯混合模型的空压机监控诊断系统及方法，提高空压机控制系统的自动化程度，实现空压机的远、近程监控和全自动监控，实现无人值守，以及实现空压机故障在线监测和诊断。

本发明是以如下技术方案实现的：一种自适应核高斯混合模型的空压机监控诊断系统，包括现场设备层、设备控制层和管理监控层，

所述的现场设备层采用PLC200作为分站，PLC200连接有用于采集水泵电机的电压、电流、运行时间，拖动电机的电压、电流、定子轴温、转子前轴温、转子后轴温、运行时间，风包内的压力、温度以及空压机的进水压力、出水压力、进水温度、出水温度的传感器，PLC20连接有用于控制空压机和水泵启和\或闭执行器；

所述的设备控制层采用PLC300作为主站，PLC300连接有上位机；主站和分站之间选用PROFIBUS-DP主从网络进行通讯，组成分布式I/O系统；主站采用PLC300为中心控制单元，基于支持向量机的变结构自适应PID控制器对整个空压机系统进行控制，上位机通过PROFIBUS-DP协议实现与主站的通讯，将每个分站的信息送到上位机，实现对整个空压机系统的监控；

所述的管理监控层过工业以太网和上位机进行通信，实现上位机的远程监控和故障诊断。

其进一步是：所述的基于支持向量机的变结构自适应PID控制器的具体算法如下：

(1)使用SVM对非线性控制对象逼近，并瞬时线性化非线性控制对象逼近模型；

(2)根据控制误差自适应调整PID控制器结构；

(3)使用多步预测思想自适应整定PID参数；

(4)根据反馈的实际压力信号与设定压力信号的误差，按照自适应PID控制算法计算控制量，实现对空压机系统的控制。

所述的空压机系统包括四台空压机，四台空压机采用主-辅-备-检修的组合方式运行，主机设备始终处于闭环控制状态。

一种自适应核高斯混合模型的空压机监控诊断方法，采用上述任一监控诊断系统，具体步骤如下：

(1)确定故障诊断所用的关键变量，采集空压机样本数据作为训练样本；

所用的监测参数包括：水泵电机的电压、电流、运行时间；拖动电机的电压、电流、定子轴温、转子前轴温、转子后轴温、运行时间；风包内的压力、温度；空压机的进水压力、出水压力、进水温度、出水温度；

(2)对采集到采集空压机样本数据进行归一化处理，将数据去掉量纲的影响；

(3)在低维子空间内使用GMM对数据集聚类揭示数据多模态；

(4)然后根据数据多模态信息，在原始数据空间利用基于熵的核函数参数估计算法寻找局部最优核函数参数；

(5)利用局部最优核函数参数使用KPCA把非线性高维多模态数据集变换低维鉴别子空间，如果算法没有达到终止条件，则转到第一步继续迭代执行，通过自适应核高斯混合模型完成对空压机系统的故障诊断。

本发明的有益效果是：提高空压机控制系统的自动化程度，实现空压机的远、近程监控和全自动监控，实现无人值守，以及实现空压机故障在线监测和诊断。

附图说明

图1是本发明的基于PLC的集控系统结构图；

图2是本发明的基于SVM的变结构自适应PID控制系统原理图；

图3是本发明压风机故障诊断系统的重要特征参数来源示意图；

图4是本发明的一种自适应核高斯混合模型的空压机故障诊断方法流程图。

具体实施方式

如图1和图2所示，一种自适应核高斯混合模型的空压机监控诊断系统，包括现场设备层、设备控制层和管理监控层，所述的现场设备层由PLC200，传感器，空压机，执行器和水泵构成，以PLC200作为分站，完成对现场设备层的控制；所述的设备控制层包括上位机和PLC300，以PLC300作为主站，主站和5个从站之间选用PROFIBUS-DP主从网络进行通讯，组成分布式I/O系统，从站的输入信号能快速地传送至主站，主站发出的指令和输出结果也能及时地送达从站执行输出，采用PLC为中心控制单元，基于支持向量机(SVM)的变结构自适应PID控制器对整个空压机系统进行控制，上位机采用西门子WinCC软件，该软件通过PROFIBUS-DP协议实现与主站的通讯，将每个分站的信息送到上位机，实现对整个空压机系统的监控；设备控制层和管理监控层通过工业以太网进行通信，从而实现上位机3的远程监控和故障诊断，而系统故障诊断的方法采用自适应核高斯混合模型。

空压机系统以PLC300为中心控制单元，使用SVM在线辨识控制模型，自适应调节PID参数，保证PID控制算法工作于最佳状态，自适应地调节空压机工作负荷，达到恒压供气的目的，提高系统的稳定性和鲁棒性；空压机通过变频软启动实现启动、停止、加载与卸荷状态转换；四台空压机采用主-辅-备-检修的组合方式运行，主机设备始终处于闭环控制状态，根据反馈的实际压力信号与设定压力信号的误差，按照自适应PID控制算法计算控制量，根据监测到的压缩空气管路压力值，自主决策控制现场空压机的开停，控制系统能够自动实现对主管路的压力控制，实现无人值守运行。

基于支持向量机的变结构自适应PID控制器的算法原理如下：

(1)使用SVM对非线性控制对象逼近，并瞬时线性化非线性对象逼近模型；

(2)根据控制误差自适应调整PID控制器结构；

(3)使用多步预测思想自适应整定PID参数；

(4)根据反馈的实际压力信号与设定压力信号的误差，按照自适应PID控制算法计算控制量，进而实现对空压机系统的控制。

上位机3采用采用西门子WinCC软件，对5个分站实时监测运行状态，同时提供故障指示及报警；界面中可以实现空压机的远程控制、运行模式切换、故障查看及空压机和水泵运行状态的查询；同时可在监控界面进行报警参数的设置，提供实时数据查询、实时报表、报警查询等的链接，可快速切换查询报表、实时数据、报警等个功能，更直观的显示出故障的所在；而且对系统运行中的电源故障、空压机脱网、停电、总送气管压力过高或过低、冷却水压过低、管路超温等故障实时监控，及时发出报警提示。

如图3，基于自适应核高斯混合模型的故障诊断方法，采用图1中的监控诊断系统，具体步骤如下：

(1)确定故障诊断所用的关键变量，采集空压机样本数据作为训练样本；

所用的监测参数包括：水泵电机的电压、电流、运行时间；拖动电机的电压、电流、定子轴温、转子前轴温、转子后轴温、运行时间；风包内的压力、温度；空压机的进水压力、出水压力、进水温度、出水温度；

(2)对采集到采集空压机样本数据进行归一化处理，具体步骤如下：

21)在低维子空间内使用GMM对数据集聚类揭示数据多模态；

22)然后根据数据多模态信息，在原始数据空间利用基于熵的核函数参数估计算法寻找局部最优核函数参数；

23)利用局部最优核函数参数使用KPCA把非线性高维多模态数据集变换低维鉴别子空间，如果算法没有达到终止条件，则转到第一步继续迭代执行，通过基于熵的自适应核高斯混合模型完成对空压机系统的故障诊断，具有重要的科学研究价值与工程应用价值。

自适应核高斯混合模型的目的就是在低维空间内对低维约简数据集进行聚类。对于高维稀疏数据集而言，GMM聚类算法往往受到成分协方差矩阵奇异的问题。如何寻找低维子空间提高数据模态之间的分开程度以及模态内部的紧凑性，是决定GMM算法性能的关键。KPCA只是保留原始高维数据集的方差信息，低维子空间只能反映的数据总体几何结构信息，无法揭示数据集隐含的几何鉴别信息。而充分利用数据集隐含鉴别信息是提高GMM性能的关键。如何选择合适的KPCA核函数参数是提高KPCA抽取数据集隐含的模态鉴别信息能力的关键。自适应核高斯混合模型通过迭代执行KPCA对高维数据维数约简以及GMM聚类算法。其中利用GMM聚类结果得到的数据鉴别信息得到KPCA最优核函数参数，从而在低维约简子空间内保持了数据的鉴别几何信息，反过来也提高GMM在低维空间的数据聚类性能。这样，通过上述迭代优化过程最终求取最优的核函数参数值以及聚类效果。如图4所示，采用熵的自适应核高斯混合模型对空压机运行状态进行故障诊断和预测，步骤如下：

1、高斯混合模型及计算过程：

高斯混合模型可以由多个单高斯分量的加权和来表示为

$p\left(x\right|\mathrm{\Θ})=\underset{i=1}{\overset{N_k}{\mathrm{\Σ}}}{\mathrm{\π}}_{i}p\left(x\right|{\mathrm{\θ}}_i)---\left(1\right)$

其中，π_i是混合参数并且满足θ_i为第i个高斯成分的参数，令X＝[x₁，x₂，…，x_N]为训练数据集合空压机运行特征的参数训练样本，假设x_i从混合分布中独立采样，我们的目标是寻找Θ使得似然函数p(X|Θ)最大。这样混合模型对应的优化问题可以写为

$L\left(\mathrm{\Θ}\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\log \underset{j=1}{\overset{N_k}{\mathrm{\Σ}}}{\mathrm{\π}}_j{p}_j\left(x\right|{\mathrm{\θ}}_j)---\left(2\right)$

$s.t.\underset{i=1}{\overset{N_k}{\mathrm{\Σ}}}{\mathrm{\π}}_i=1$

上述优化是非凸优化问题，为了简化上述似然表示形式，引入隐含指示变量z_ij∈{0，1}，且z_ij＝1。表示x_i采样于分布分量j，反之则z_ij＝0；

那么p(Z|Θ)表示为

$p\left(Z\right|\mathrm{\Θ})=\underset{i=1}{\overset{N}{\mathrm{\Π}}}p\left(z_i\right|\mathrm{\Θ})=\underset{i=1}{\overset{N}{\mathrm{\Π}}}\underset{j=1}{\overset{N_k}{\mathrm{\Π}}}{{\mathrm{\π}}_j}^{z_{\mathrm{ij}}}---\left(3\right)$

关于Θ和的条件分布的似然函数变为

$p\left(X\right|\mathrm{\Θ},Z)=\underset{i=1}{\overset{N}{\mathrm{\Π}}}\underset{j=1}{\overset{N_k}{\mathrm{\Π}}}p{\left(x_i\right|{\mathrm{\θ}}_j)}^{z_{\mathrm{ij}}}---\left(4\right)$

因此，对数完整似然函数变为

$L\left(\mathrm{\Θ}\right)=\log p(X,Z|\mathrm{\Θ})=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_k}{\mathrm{\Σ}}}{z}_{\mathrm{ij}}\log {\mathrm{\π}}_{j}p\left(x_i\right|{\mathrm{\θ}}_j)---\left(5\right)$

2、本发明提出的自适应求解最优参数N_i的过程：

最优参数可以通过期望最大化算法等方法求取。高阶统计量skewness和kurtosis能够度量高斯混合密度拟合特定类的训练数据[]。令N_i为第i成分的训练数据数量，和分别为第i个成分关于第d维的统计量skewness和kurtosis，其计算如下

${\hat{s}}_{i,d}=\frac{\frac{1}{N_i}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{h}_i^n{(x_{\mathrm{id}}^n-{\widehat{\mathrm{\μ}}}_{\mathrm{id}})}^3}{{\widehat{\mathrm{\σ}}}_{\mathrm{id}}^3}---\left(6\right)$

${\hat{k}}_{i,d}=\frac{\frac{1}{N_i}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{h}_i^n{(x_{\mathrm{id}}^n-{\widehat{\mathrm{\μ}}}_{\mathrm{id}})}^4}{{\widehat{\mathrm{\σ}}}_{\mathrm{id}}^4}---\left(7\right)$

其中，是属于第i成分的后验概率，和第i个成分关于第d维的均值和标准差。当第i个成分关于第d维的分布是高斯分布，Skewness和kurtosis估计逼近0，否则高阶统计量偏离0。而关于成分i的skewness和kurtosis分布可以计算为

${\hat{s}}_i=\frac{1}{D}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}\left|{\hat{s}}_{i,d}\right|---\left(8\right)$

${\hat{k}}_i=\frac{1}{D}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}\left|{\hat{k}}_{i,d}\right|---\left(9\right)$

关于GMM的非高斯性可以定义为

$\mathrm{\Φ}=\underset{c=1}{\overset{C}{\mathrm{\Σ}}}{\mathrm{\π}}_c({\hat{s}}_c+{\hat{k}}_c)---\left(10\right)$

该度量估计GMM对训练数据本质分布拟合的程度，使用下面准则选择拟合度最差的成分

$c^{*}=\underset{c}{\arg \max }{\mathrm{\π}}_c({\hat{s}}_c+{\hat{k}}_c)---\left(11\right)$

假设第c个成分非高斯性最大，使用文[]非负对该成分进行分裂，新的成分参数初始化为Φ＝{Θ_k1，Θ_k2}。分裂过程为[]

$a_{c_1}=u_1{a}_{c^{*}}---\left(12\right)$

$a_{c_2}=(1-u_1)a_{c^{*}}---\left(13\right)$

${\mathrm{\μ}}_{c_1}={\mathrm{\μ}}_{c^{*}}-\left(\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{u}_2^i\sqrt{{\mathrm{\λ}}_{c^{*}}^i}{v}_{c^{*}}^i\right)\sqrt{\frac{a_{c_1}}{a_{c_2}}}---\left(14\right)$

${\mathrm{\μ}}_{c_2}={\mathrm{\μ}}_{c^{*}}+\left(\underset{i=1}{\overset{d}{\mathrm{\Σ}}}{u}_2^i\sqrt{{\mathrm{\λ}}_{c^{*}}^i}{v}_{c^{*}}^i\right)\sqrt{\frac{a_{c_1}}{a_{c_2}}}---\left(15\right)$

${\mathrm{\Λ}}_1=\mathrm{diag}\left(u3\right)\mathrm{diag}(\mathrm{\ι}-u2)\mathrm{diag}(\mathrm{\ι}+u2){\mathrm{\Λ}}_{*}\frac{a_{c^{*}}}{a_{c_1}}---\left(16\right)$

${\mathrm{\Λ}}_2=\mathrm{diag}(\mathrm{\ι}-u3)\mathrm{diag}(\mathrm{\ι}-u2)\mathrm{diag}(\mathrm{\ι}+u2){\mathrm{\Λ}}_{*}\frac{a_{c^{*}}}{a_{c_2}}---\left(17\right)$

$V_1=D{V}_{c^{*}}---\left(18\right)$

$V_2=D^T{V}_{c^{*}}---\left(19\right)$

这里和分别是矩阵的升序特征值和对应的特征向量。是对角矩阵，D表示D×D旋转矩阵，矩阵的每一列为正交单位向量。该矩阵由D(D-1)/2个不同符合均匀分布U(0，1)的数据独立地产生下三角矩阵建立的。ι是元素全为1的D×1向量，和表示2D+1随机变量，它们由下面方法产生

$u_1~\mathrm{Beta}\left(\mathrm{2,2}\right),u_2^1~\mathrm{Beta}\left(\mathrm{1,2}D\right),u_3^1~\mathrm{Beta}(1,D)---\left(20\right)$

$u_2^j~U(-\mathrm{1,1})\mathrm{and}{u}_3^j~U\left(\mathrm{0,1}\right)\mathrm{with\; j}=\mathrm{2,3},...,D---\left(21\right)$

这里Beta(·，·)是beta分布。

3、本发明提出的求解最优参数σ和w_d参数过程：

核主元分析(KPCA)的最大的优势是允许我们使用一个mercer核函数计算高维空间上的点积。因此，无需从输入空间的数据投影到高维空间的映射形式，并且有效地描述了数据的非线性结构。KPCA方法是一种多任务的有效预处理方法，因其物理意义简单、便于应用等优点引起广泛关注。

具体地，给定N个训练样本x₁，x₂，…，x_N，令φ：x→φ(x)∈H为从输入空间到高维特征空间H的非线性映射，使用与重建核希尔伯特空间(Reproducing Kernel Hilbert Space，RKHS)相关的核函数计算特征空间H中内积，其表示为k(x，x′)＝φ(x)^Tφ(x′)。两种常用的核函数为多项式核函数和Gaussian径向基核函数，其定义分别为

$k(x_i,x_j)={(x_j^T{x}_i+1)}^d---\left(22\right)$

$k(x_i,x_j)=\exp (-\frac{{\left|\right|x_j-x_i\left|\right|}^2}{\mathrm{\γ}}),\mathrm{\γ}\∈R^{+}---\left(23\right)$

特征空间的映射样本的均值由下式计算为

$\stackrel{\‾}{\mathrm{\φ}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{\φ}\left(x_i\right)=\mathrm{\ΦS}---\left(24\right)$

这里Φ＝[φ(x₁)，φ(x₂)，…，φ(x_N)]，S＝1_N/N，1_N表示全部元素为1的N×1向量。特征空间中样本的中心化方法为

$\widehat{\mathrm{\φ}}\left(x_i\right)=\mathrm{\φ}\left(x_1\right)-\stackrel{\‾}{\mathrm{\φ}},i=\mathrm{1,2},...,N---\left(25\right)$

那么，由下式构建特征空间的协方差矩阵的中心化核矩阵为

$\stackrel{\‾}{K}=K-\frac{1}{N}{1}_{\mathrm{NN}}K-\frac{1}{N}K{1}_{\mathrm{NN}}+\frac{1}{N^2}{1}_{\mathrm{NN}}K{1}_{\mathrm{NN}}---\left(26\right)$

其中，1_NN表示N×N的单位矩阵，中心化核矩阵函数的第(i，j)个元素为 $\stackrel{\‾}{k}(x_i,x_j)=\stackrel{\‾}{k}(x_i,x_j)-\frac{1}{N}{1}_{1N}{k}_{x_i}-\frac{1}{N}{1}_{1N}{k}_{x_j}+\frac{1}{N^2}{1}_{1N}K{1}_{N1},$ $k_{x_i}={[k(x_i,x_1),k(x_i,x_2),...,k(x_i,x_N)]}^T$ 和K＝[k(x_i，x_j)]_{i，j＝1，2，…，N}。由于核函数是对称、连续和正定的，那么矩阵K也是对称和半正定的。KPCA由下式表示的特征值问题求取，即

$\stackrel{\‾}{K}{\mathrm{\α}}_i={\mathrm{\λ}}_i{\mathrm{\α}}_i---\left(27\right)$

这里λ_i是的第i个最大特征值，α_i为对应λ_i的特征向量。由于每个成分是特征空间的协方差矩阵的特征向量，每个成分由训练样本张成，其相应的系数为α的坐标。因此映射样本投影到第i个成分v_i由下式给出

${\mathrm{\β}}_i=\stackrel{\‾}{\mathrm{\φ}}{\left(x\right)}^T{v}_i=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\α}}_{\mathrm{in}}k(x,x_n)={\mathrm{\α}}_i^T{\stackrel{\‾}{k}}_x---\left(28\right)$

假设，使用矩阵的d个最大特征值对应的d个成分进行维数约简，到d个成分后的投影表示为

$P_d\left(x\right)={[{\mathrm{\α}}_1,{\mathrm{\α}}_2,...,{\mathrm{\α}}_d]}^T{\stackrel{\‾}{k}}_x=A^T{\stackrel{\‾}{k}}_x---\left(29\right)$

这里A＝[α₁，α₂，…，α_d].在实际应用中，KPCA的优越性能取决于核函数参数。常规的KPCA模型参数即尺度参数只是一个正数的作用是控制核映射的非线性，并且均等对待各个变量。然而不同变量对揭示低维子空间中数据结构的作用不同，因此只使用单个尺度参数往往是不够的。为此，该发明提出的改进型高斯核函数为

$k(x,y)=\exp \left(\frac{{(x-y)}^{T}W(x-y)}{2}\right)=\exp (-\frac{1}{2{\mathrm{\σ}}^2}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_d{(x_d-y_d)}^2)---\left(30\right)$

对大多数监督和半监督学习方法，使用必须连接(ML)和不可能连接(CL)。约束表示的监督信息可以有效提高学习算法的性能。其中，ML约束表示两个观测样本必须属于相同聚类或类别，CL约束指定不属于同类的两个观测值。

令φ：x→φ(x)∈H为从输入空间到高维特征空间H的非线性映射几何，空间H的点积由定义为k(x，x′)＝φ(x)^Tφ(x′)的Mercer核函数计算。常用的高斯核函数定义为本发明采用的高斯核函数定义为

$k(x,y)=\exp \left(\frac{{(x-y)}^{T}W(x-y)}{2}\right)=\exp (-\frac{1}{2{\mathrm{\σ}}^2}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_d{(x_d-y_d)}^2)---\left(31\right)$

这里W＝σ^-2diag(w₁，w₂，…，w_D).w_d为数据第d维的权重。从概率论的观点看，w_d可以看作第d维数据对数据集贡献的概率，因而可理解为维度权重。

给出ML和CL，最优核函数参数可以通过使用CL和ML求取。优化问题可以定义为

$\min L(W,{\mathrm{\σ}}^2)=-\underset{(x_i,x_j)\∈\mathrm{ML}}{\mathrm{\Σ}}{\left|\right|\mathrm{\φ}\left(x_i\right)-\mathrm{\φ}\left(x_j\right)\left|\right|}^2+\underset{(x_i,x_j)\∈\mathrm{CL}}{\mathrm{\Σ}}{\left|\right|\mathrm{\φ}\left(x_i\right)-\mathrm{\φ}\left(x_j\right)\left|\right|}^2+\mathrm{\η}\underset{i=1}{\overset{D}{\mathrm{\Σ}}}{w}_i\ln w_i---\left(32\right)$

$s.t.\underset{i=1}{\overset{D}{\mathrm{\Σ}}}{w}_i=\mathrm{1,0}\≤w_i\≤1---\left(33\right)$

最小化等价于最大化一个聚类和其他聚类之间的聚类间相似度，同时最小化聚类的类内相似度.最后项使用目标函数的关于维度权重作为正则化项，目的用于避免聚类过程的过拟合问题。正参数γ控制类内相似度和类间相似度之间的折中。目标函数中正则化项的目的是使得更多的维度辨识聚类，解决稀疏数据聚类只使用少量维度聚类的问题。因此在存在冗余特征时，在优化问题显式刻画稀疏性，较好地提高KPCA性能。所在该策略较好解决聚类高维、稀疏数据问题。

使用核技巧并丢弃不影响最终优化结果的常数项，上述优化问题可以转化文如下优化问题

$\min J(W,{\mathrm{\σ}}^2)=-\underset{(x_i,x_j)\∈\mathrm{CL}}{\mathrm{\Σ}}k(x_i,x_j)-\underset{(x_i,x_j)\∈\mathrm{ML}}{\mathrm{\Σ}}k(x_i,x_j)+\mathrm{\η}\underset{i=1}{\overset{D}{\mathrm{\Σ}}}{w}_i\ln w_i---\left(34\right)$

$s.t.\underset{i=1}{\overset{D}{\mathrm{\Σ}}}{w}_i=1---\left(35\right)$

关于w_i和σ²的最优解不能通过直接求解上述优化问题求得。这里采用已有的多乘子方法求解上述优化问题。上述优化问题对应的Lagrange乘子函数可以写为

$L(W,{\mathrm{\σ}}^2,\mathrm{\λ})=J(W,{\mathrm{\σ}}^2)-\mathrm{\λ}(\underset{i=1}{\overset{D}{\mathrm{\Σ}}}{w}_i-1)---\left(36\right)$

广义目标函数W，σ²，δ和λ^*分别可以在下式表达为

$\arg \min M(W,{\mathrm{\σ}}^2,\mathrm{\δ},\mathrm{\λ})=L(W,{\mathrm{\σ}}^2,\mathrm{\λ})+\frac{\mathrm{\δ}}{2}{(\underset{i=1}{\overset{D}{\mathrm{\Σ}}}{w}_i-1)}^2---\left(37\right)$

为了计算上述无约束优化问题，在固定λ和δ的条件下，用梯度下降算法计算W和σ²，计算过程如下：

$\frac{\∂M}{\∂\mathrm{\σ}}=\frac{\∂J}{\∂\mathrm{\σ}}=\underset{(x_i,x_j)\∈\mathrm{CL}}{\mathrm{\Σ}}\frac{\∂k(x_i,x_j)}{\∂\mathrm{\σ}}-\underset{(x_i,x_j)\∈\mathrm{ML}}{\mathrm{\Σ}}\frac{\∂k(x_i,x_j)}{\∂\mathrm{\σ}}---\left(38\right)$

$\frac{\∂M}{\∂w_i}=\frac{\∂L}{\∂w_i}+\mathrm{\δ}(\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_d-1)---\left(39\right)$

式中，

$\frac{\∂k(x_i,x_j)}{\∂\mathrm{\σ}}=\left(\frac{1}{{\mathrm{\σ}}^3}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_d{(x_d-y_d)}^2\right)\exp (-\frac{1}{2{\mathrm{\σ}}^2}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_d{(x_d-y_d)}^2)---\left(40\right)$

$\frac{\∂L}{\∂w_i}=\underset{(x_i,x_j)\∈\mathrm{CL}}{\mathrm{\Σ}}\frac{\∂k(x_i,x_j)}{\∂w_i}-\underset{(x_i,x_j)\∈\mathrm{ML}}{\mathrm{\Σ}}\frac{\∂k(x_i,x_j)}{\∂w_i}+1+\ln w_i-\mathrm{\λ}---\left(38\right)$

$\frac{\∂k(x_i,x_j)}{\∂w_i}=-\frac{1}{2{\mathrm{\σ}}^2}{(x_i-y_i)}^2\exp (-\frac{1}{2{\mathrm{\σ}}^2}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_d{(x_d-y_d)}^2)---\left(42\right)$

因此，可以得到

$\mathrm{\σ}(t+1)=\mathrm{\σ}\left(t\right)-\mathrm{\ρ}\frac{\∂M}{\∂\mathrm{\σ}}|\mathrm{\σ}=\mathrm{\σ}\left(t\right)---\left(43\right)$

$w_i(t+1)=w_i\left(t\right)-\mathrm{\ρ}\frac{\∂M}{\∂w_i}|w_i=w_i\left(t\right)---\left(44\right)$

式中，单位步长参数ρ可以通过现行搜索算法得到。利用梯度下降算法可以保证局部参数σ和w_d的最优。

$\mathrm{\λ}(t+1)=\mathrm{\λ}\left(t\right)-\mathrm{\δ}(\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_d-1)---\left(45\right)$

$\mathrm{\δ}(t+1)=\left\{\begin{array}{cc}{c}^{\′}\mathrm{\δ}\left(t\right),& \mathrm{if}|\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_d\left(t\right)-1|/|\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_d(t-1)-1|>{\mathrm{\ϵ}}^{\mathrm{\δ}}\\ \mathrm{\δ}\left(t\right),& \mathrm{otherwise}\end{array}\right.---\left(46\right)$

式中，c′＞1代表放大系数，ε^δ∈(0，1).这些参数W，σ²，λ，δ通过迭代更新直到达到终止条件乘法优化的细节可以参考[]。

4、建立压风机的检测指标

基于GMM的设备监测方法，首先根据训练数据估计EKAGMM的参数，并计算监测数据对每个高斯成分(或聚类)的后验概率。然后根据局部后验概率计算全局监测推理指标进行故障检测。样本x_i的条件概率密度函数表示为

$p\left(x_i\right|w_c,o_c,{\mathrm{\σ}}_c^2)={\left(\underset{d=1}{\overset{D}{\mathrm{\Π}}}\frac{w_{\mathrm{cd}}}{2\mathrm{\π}{\mathrm{\σ}}_c^2}\right)}^{-\frac{1}{2}}\exp (-\frac{1}{2{\mathrm{\σ}}_c^2}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_{\mathrm{td}}\left({(x_{\mathrm{id}}-o_{\mathrm{cd}})}^2\right))---\left(47\right)$

根据贝叶斯准则，样本x_i的后验概率可以通过下式计算

$p\left(c\right|x)=\frac{a_c\left(\underset{d=1}{\overset{D}{\mathrm{\Π}}}\sqrt{\frac{w_{\mathrm{cd}}}{{\mathrm{\σ}}_c^2}}\right)\exp (-\frac{1}{2{\mathrm{\σ}}_c^2}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_{\mathrm{cd}}\left({(x_{\mathrm{id}}-o_{\mathrm{cd}})}^2\right))}{\underset{t=1}{\overset{C}{\mathrm{\Σ}}}\left(\underset{d=1}{\overset{D}{\mathrm{\Π}}}\sqrt{\frac{w_{\mathrm{cd}}}{{\mathrm{\σ}}_t^2}}\right)\exp (-\frac{1}{2{\mathrm{\σ}}_t^2}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_{\mathrm{td}})\left({(x_{\mathrm{id}}-o_{\mathrm{cd}})}^2\right)}---\left(48\right)$

假设每个高斯成分服从单一模态高斯分布，那么从一个高斯成分采样得到的样本x_t的马氏距离(MD)服从分布。表示样本x_i与高斯成分MD平方，这里自由度D等于输入样本为维数， $D_{\mathrm{local}}^c\left(x_i\right)={(x_i-o_c)}^T{\mathrm{\Σ}}_c^{-1}(x_i-o_c)=\frac{1}{{\mathrm{\σ}}_c^2}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_{\mathrm{cd}}\left({(x_{\mathrm{id}}-o_{\mathrm{cd}})}^2\right).$ 基于局部MD的监测样本属于每个聚类的概率可以由下式估计

$p^{\left(c\right)}\left(x_t\right)=\mathrm{Pr}\{\left(D_{\mathrm{local}}^c\left(x\right)\right|x\∈{\mathrm{\Θ}}_c)\≤\left(D_{\mathrm{local}}^c\left(x_t\right)\right|x_t\∈{\mathrm{\Θ}}_c)\}---\left(49\right)$

p^(c)(x_t)也可以由带有合适自由度的χ²概率密度函数计算，它表示样本x_t属于成分Θ_c的概率。类似于基于Bayesian推理的推理监测统计量构建方法，基于EAKGMM的推理监测统计量(EAKGMMIMS)由下式建立

$\begin{array}{c}\mathrm{EAKGMMIMS}\left(x_t\right)=\underset{c=1}{\overset{C}{\mathrm{\Σ}}}{p}^{\left(c\right)}\left(x_t\right)p\left(c\right|x_t)\\ =\underset{c=1}{\overset{C}{\mathrm{\Σ}}}{p}^{\left(c\right)}\left(x_t\right)\frac{a_c\left(\underset{d=1}{\overset{D}{\mathrm{\Π}}}\sqrt{\frac{w_{\mathrm{cd}}}{{\mathrm{\σ}}_c^2}}\right)\exp (-\frac{1}{2{\mathrm{\σ}}_c^2}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_{\mathrm{cd}}\left({(x_{\mathrm{id}}-o_{\mathrm{cd}})}^2\right))}{\underset{t=1}{\overset{C}{\mathrm{\Σ}}}\left(\underset{d=1}{\overset{D}{\mathrm{\Π}}}\sqrt{\frac{w_{\mathrm{cd}}}{{\mathrm{\σ}}_t^2}}\right)\exp (-\frac{1}{2{\mathrm{\σ}}_t^2}\underset{d=1}{\overset{D}{\mathrm{\Σ}}}{w}_{\mathrm{td}}\left({(x_{\mathrm{id}}-o_{\mathrm{cd}})}^2\right))}\end{array}---\left(50\right)$

全局监测统计量利用所有高斯成分的后验概率，表示观测样本在正常状态内的概率。由于无需确定监测样本x_t属于的单个高斯成分，全局监测统计量指标能够避免由误分类引起的误监测的潜在风险。由于p^(c)(x_t)和p(c|x_t)均小于1，很容易可知EAKGMMIMS(x_t)也小于1.因此可以通过EAKGMMIMS是否超出(1-a)100％判定设备状态是否正常还是故障状态。

6、故障诊断和分类

该方法可以通过比较各个监测变量的贡献率来确定其故障源。在第c类中第d个被监测变量被估计通过一阶偏导数

$C_{D_c}^{\left(d\right)}\left(x\right)=\frac{\∂D_{\mathrm{local}}^c\left(x\right)}{\∂x_d}=\frac{2}{{\mathrm{\σ}}_c^2}{w}_{\mathrm{cd}}\left((x_{\mathrm{id}}-o_{\mathrm{cd}})\right)---\left(51\right)$

然后，EEGMMIMS的d维监测变量的总贡献率可以通过下式得到

$C_{\mathrm{EAKGMMIMS}}^{\left(d\right)}\left(x\right)=\underset{c=1}{\overset{C}{\mathrm{\Σ}}}{p}^{\left(c\right)}\left(x\right)C_{D_c}^{\left(d\right)}\left(x\right)---\left(52\right)$

控制上限的显著性水平a可以通过下式得到

$\mathrm{DL}=\frac{D(N-1)}{N-D}{F}_{D,N-D,a}---\left(53\right)$

因此，度量服从近似χ²分布。

本发明能够较好地跟踪设定的压力值，解决了空压机负荷大范围变化的问题，提高了控制精度，具有较好的鲁棒性能；克服了现有空压机监测算法没有很好考虑变量之间的相关性而导致误报率过高的问题，并且适应设备多工况运行状态监测。采用了基于自适应核高斯混合模型对运行故障进行实时的监测与诊断，具有更高的诊断精度。

Claims (4)

1.一种自适应核高斯混合模型的空压机监控诊断系统，包括现场设备层、设备控制层和管理监控层，其特征在于：

所述的现场设备层采用PLC200作为分站，PLC200连接有用于采集水泵电机的电压、电流、运行时间，拖动电机的电压、电流、定子轴温、转子前轴温、转子后轴温、运行时间，风包内的压力、温度以及空压机的进水压力、出水压力、进水温度、出水温度的传感器，PLC200连接有用于控制空压机和水泵启和\或闭执行器；

所述的设备控制层采用PLC300作为主站，PLC300连接有上位机；主站和分站之间选用PROFIBUS-DP主从网络进行通讯，组成分布式I/O系统；主站采用PLC300为中心控制单元，基于支持向量机的变结构自适应PID控制器对整个空压机系统进行控制，上位机通过PROFIBUS－DP协议实现与主站的通讯，将每个分站的信息送到上位机，实现对整个空压机系统的监控；

所述的管理监控层通过工业以太网和上位机进行通信，实现上位机的远程监控和故障诊断。

2.根据权利要球1所述的一种自适应核高斯混合模型的空压机监控诊断系统，其特征在于：所述的基于支持向量机的变结构自适应PID控制器的具体算法如下：

(1)使用SVM对非线性控制对象逼近，并瞬时线性化非线性控制对象逼近模型；

(2)根据控制误差自适应调整PID控制器结构；

(3)使用多步预测思想自适应整定PID参数；

(4)根据反馈的实际压力信号与设定压力信号的误差，按照自适应PID控制算法计算控制量，实现对空压机系统的控制。

3.根据权利要球1或2所述的一种自适应核高斯混合模型的空压机监控诊断系统，其特征在于：所述的空压机系统包括四台空压机，四台空压机采用主-辅-备-检修的组合方式运行，主机设备始终处于闭环控制状态。

4.一种自适应核高斯混合模型的空压机监控诊断方法，其特征在于：采用1-3中任一监控诊断系统，具体步骤如下：

(1)确定故障诊断所用的关键变量，采集空压机样本数据作为训练样本；

所用的监测参数包括：水泵电机的电压、电流、运行时间；拖动电机的电压、电流、定子轴温、转子前轴温、转子后轴温、运行时间；风包内的压力、温度；空压机的进水压力、出水压力、进水温度、出水温度；

(2)对采集到采集空压机样本数据进行归一化处理，将数据去掉量纲的影响；

(3)在低维子空间内使用GMM对数据集聚类揭示数据多模态；

(4)然后根据数据多模态信息，在原始数据空间利用基于熵的核函数参数估计算法寻找局部最优核函数参数；

(5)利用局部最优核函数参数使用KPCA把非线性高维多模态数据集变换低维鉴别子空间，如果算法没有达到终止条件，则转到第一步继续迭代执行，通过自适应核高斯混合模型完成对空压机系统的故障诊断。