CN110083860B - 一种基于相关变量选择的工业故障诊断方法 - Google Patents

一种基于相关变量选择的工业故障诊断方法 Download PDF

Info

Publication number
CN110083860B
CN110083860B CN201910189841.4A CN201910189841A CN110083860B CN 110083860 B CN110083860 B CN 110083860B CN 201910189841 A CN201910189841 A CN 201910189841A CN 110083860 B CN110083860 B CN 110083860B
Authority
CN
China
Prior art keywords
fault
spe
industrial
matrix
vector
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910189841.4A
Other languages
English (en)
Other versions
CN110083860A (zh
Inventor
邓瑞祥
张颖伟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northeastern University China
Original Assignee
Northeastern University China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northeastern University China filed Critical Northeastern University China
Priority to CN201910189841.4A priority Critical patent/CN110083860B/zh
Publication of CN110083860A publication Critical patent/CN110083860A/zh
Application granted granted Critical
Publication of CN110083860B publication Critical patent/CN110083860B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation

Abstract

本发明涉及故障监测与诊断技术领域，提供一种基于相关变量选择的工业故障诊断方法，首先采集正常工况样本和每种工业故障的故障工况样本，然后根据正常工况样本建立PCA监测模型，再基于LASSO算法，同时考虑正常工况数据与故障工况数据以及从正常工况数据到故障工况数据的相对变化来建立每种工业故障的故障相关变量选择模型，通过结合故障重构思想的修正的LARS算法求解模型以精准地选择出每种故障的相关变量，接着基于故障子空间的思想建立每个子空间的PCA监测模型，最后实时采集工业生产样本并对工业故障进行在线实时检测与诊断。本发明能够在检测出工业故障的同时对故障类型进行诊断，提高工业故障诊断的效率及准确性。

Description

一种基于相关变量选择的工业故障诊断方法
技术领域
本发明涉及故障监测与诊断技术领域，特别是涉及一种基于相关变量选择的工业故障诊断方法。
背景技术
随着工业过程变得更加集成化、复杂化和智能化，如何有效地对生产过程进行监控从而对故障进行检测并诊断，成为保证生产安全、提升产品质量及节约能源的关键。随着传感器、自动化和计算机技术的大规模应用，工业过程中可测得的过程变量的个数快速增长，极大地促进了基于数据的工业过程监测方法的发展。其中，使用最多的工业过程监测方法是多变量统计过程监测(Multivariate Statistical Process Monitoring,MSPM)方法。作为应用最多的MSPM方法，主成分分析(Principal Component Analysis，PCA)法利用矩阵变换的思想从高度相关的高维原始数据中提取相互独立的低维主元，然后通过在低维子空间中计算相应的统计量完成工业过程监测。
然而，工业过程中可测得变量个数的增多，为检测到故障后的后续处理工作带来了困难；因为一个典型的工业生产过程往往包含大量的生产环节和回路，故障往往局限于其中的一小部分，在检测到故障的发生后，要对故障发生的位置及类型进行快速诊断就需要耗费相关领域专家大量的时间与精力。现有的方法不能在检测出故障的同时提取故障相关变量进行根源分析并诊断故障类型，从而故障诊断的效率较低且准确性不高。
发明内容
针对现有技术存在的问题，本发明提供一种基于相关变量选择的工业故障诊断方法，能够精准地提取每种工业故障的相关变量以进行故障根源分析，能够在检测出工业故障的同时对故障类型进行诊断，提高工业故障诊断的效率及准确性。
本发明的技术方案为：
一种基于相关变量选择的工业故障诊断方法，其特征在于，包括下述步骤：
步骤1：在工业生产过程中，确定m个可测得变量{a1,a2,…,am}，对正常工况下n个时刻的m个可测得变量的数据进行采集，得到n个正常工况样本，构成初始正常工况矩阵，对初始正常工况矩阵进行标准化处理，得到正常工况矩阵X＝(aij)n×m，aij为第i个正常工况样本中第j个可测得变量的值，i∈{1,2,…,n}，j∈{1,2,…,m}；确定F种工业故障，对每种工业故障下n个时刻的m个可测得变量的数据进行采集，得到每种工业故障的n个故障工况样本，构成每种工业故障的初始故障工况矩阵，对每种工业故障的初始故障工况矩阵进行标准化处理，得到第f种工业故障的故障工况矩阵Xf＝(aij f)n×m，f∈{1,2,...,F}，aij f为第f种工业故障的第i个故障工况样本中第j个可测得变量的值；
步骤2：利用正常工况矩阵X建立PCA监测模型为
Figure BDA0001994054830000021
其中，P∈Rm×l为负载矩阵，l为矩阵X的主元个数，I∈Rm×m为单位矩阵，
Figure BDA0001994054830000022
为正常工况矩阵X在主元空间中的投影，E为正常工况矩阵X的残差矩阵；
计算每个正常工况样本x＝(ai1,ai2,…,aim)T的T2统计量和SPE统计量分别为
Figure BDA0001994054830000028
SPE＝||(I-PPT)x||2
其中，t＝PTx∈Rl×1为得分向量，
Figure BDA0001994054830000023
将n个正常工况样本的T2统计量和SPE统计量分别存入向量T2∈Rn×1和向量SPE∈Rn ×1中，并计算正常工况样本的T2统计量的置信度为α的上限
Figure BDA0001994054830000024
和SPE统计量的置信度为α的上限SPEα；其中，α＝95％；
步骤3：通过负载矩阵P对第f种工业故障的故障工况矩阵Xf进行分解：
Xf＝XfPPT+Xf(I-PPT)
计算第f种工业故障的每个故障工况样本的T2统计量和SPE统计量，并将第f种工业故障的n个故障工况样本的T2统计量和SPE统计量分别存入向量
Figure BDA0001994054830000025
和向量SPEf∈Rn×1中；
步骤4：计算第f种工业故障的故障工况矩阵相对于正常工况矩阵的变化量矩阵为
ΔX＝Xf-X
计算第f种工业故障的向量
Figure BDA0001994054830000026
相对于向量T2的变化量矩阵为
Figure BDA0001994054830000027
计算第f种工业故障的向量SPEf相对于向量SPE的变化量矩阵为
ΔSPE＝SPEf-SPE
步骤5：基于最小绝对值收敛及选择算子回归算法LASSO，构建第f种工业故障的故障相关变量选择模型为
Figure BDA0001994054830000031
Figure BDA0001994054830000032
其中，
Figure BDA0001994054830000033
和βSPE分别为T2统计量和SPE统计量的回归系数向量，
Figure BDA0001994054830000034
和μSPE分别为T2统计量和SPE统计量的惩罚系数，
Figure BDA00019940548300000331
和μSPE分别用来控制
Figure BDA0001994054830000035
和βSPE的稀疏性；
步骤6：基于修正的最小角回归算法LARS，求解第f种工业故障的故障相关变量选择模型，具体步骤如下：
步骤6.1：对于第f种工业故障的故障相关变量选择模型的第一个子模型
Figure BDA0001994054830000036
确定第f种工业故障的在T2统计量下的估计向量
Figure BDA0001994054830000037
和第f种工业故障的在T2统计量下的故障相关变量集合
Figure BDA0001994054830000038
的初始值分别为
Figure BDA0001994054830000039
Figure BDA00019940548300000310
步骤6.2：计算相关性向量为
Figure BDA00019940548300000311
计算相关性向量
Figure BDA00019940548300000312
中元素的最大值为
Figure BDA00019940548300000313
将i添加至故障相关变量集合
Figure BDA00019940548300000314
中：
Figure BDA00019940548300000315
其中，
Figure BDA00019940548300000316
为相关性向量
Figure BDA00019940548300000317
的第i个元素，i∈{1,2,…,m}；
步骤6.3：令
Figure BDA00019940548300000318
构建矩阵
Figure BDA00019940548300000319
为
Figure BDA00019940548300000320
其中，
Figure BDA00019940548300000321
为
Figure BDA00019940548300000322
的符号，
Figure BDA00019940548300000323
为矩阵ΔX的第i列向量，
Figure BDA00019940548300000324
为全1向量，
Figure BDA00019940548300000325
为故障相关变量集合
Figure BDA00019940548300000326
中的元素个数；
步骤6.4：计算估计向量
Figure BDA00019940548300000327
的更新方向
Figure BDA00019940548300000328
和更新步长
Figure BDA00019940548300000329
分别为：
Figure BDA00019940548300000330
Figure BDA0001994054830000041
其中，
Figure BDA0001994054830000042
为故障相关变量集合
Figure BDA0001994054830000043
的补集，aj为向量
Figure BDA0001994054830000044
中的第j个元素；
步骤6.5：更新估计向量
Figure BDA0001994054830000045
步骤6.6：采用上述步骤6.1至步骤6.5的方法，对第f种工业故障的故障相关变量选择模型的第二个子模型
Figure BDA0001994054830000046
进行求解，得到第f种工业故障的在SPE统计量下的估计向量
Figure BDA00019940548300000419
和第f种工业故障的在SPE统计量下的故障相关变量集合ΓSPE；
步骤6.7：对集合
Figure BDA00019940548300000418
与ΓSPE求取并集，得到第f种工业故障的故障相关变量集合为
Figure BDA0001994054830000048
步骤6.8：利用故障相关变量集合Γf和T2统计量，对第f种工业故障的每个故障工况样本xf＝(ai1 f,ai2 f,…,aim f)T进行重构，建立第一重构模型为
Figure BDA0001994054830000049
其中，
Figure BDA00019940548300000410
为故障方向矩阵，
Figure BDA00019940548300000411
为第f种工业故障的故障相关变量集合Γf中元素的个数也即故障相关变量的个数，矩阵Ξ的第i行第j列的元素为
Figure BDA00019940548300000412
(Γf)j为集合Γf中的第j个元素；
Figure BDA00019940548300000413
为第一幅值估计向量；
求解第一重构模型，得到当前最优的第一幅值估计向量为
Figure BDA00019940548300000414
从而，得到第f种工业故障的第一重构样本为
Figure BDA00019940548300000415
步骤6.9：利用故障相关变量集合Γf和SPE统计量，对第f种工业故障的每个故障工况样本xf＝(ai1 f,ai2 f,…,aim f)T进行重构，建立第二重构模型为
Figure BDA00019940548300000416
其中，
Figure BDA00019940548300000417
为第二幅值估计向量；
求解第二重构模型，得到当前最优的第二幅值估计向量为
hSPE＝(ΞT(I-PPT)Ξ)-1ΞT(I-PPT)xf
从而，得到第f种工业故障的第二重构样本为
zSPE＝xf-ΞhSPE
＝xf-Ξ(ΞT(I-PPT)Ξ)-1ΞT(I-PPT)xf
步骤6.10：计算第f种工业故障的每个第一重构样本
Figure BDA0001994054830000051
的T2统计量为
Figure BDA0001994054830000052
其中，
Figure BDA0001994054830000053
为样本
Figure BDA0001994054830000054
的得分向量；
计算第f种工业故障的每个第二重构样本zSPE∈Rm×1的SPE统计量为
Figure BDA0001994054830000055
将第f种工业故障的n个第一重构样本的T2统计量和n个第二重构样本的SPE统计量分别存入向量
Figure BDA0001994054830000056
和向量
Figure BDA0001994054830000057
中；
步骤6.11：将向量
Figure BDA0001994054830000058
中的每个元素与上限
Figure BDA0001994054830000059
进行比较、将向量
Figure BDA00019940548300000510
中的每个元素与上限SPEα进行比较：
若向量
Figure BDA00019940548300000511
中有98％以上的元素低于上限
Figure BDA00019940548300000512
且向量
Figure BDA00019940548300000513
中有98％以上的元素低于上限SPEα，则认为重构样本的T2统计量和SPE统计量均转为正常状态，从而第f种工业故障的所有故障相关变量都已经位于集合Γf中，进入步骤7；
否则，重复上述步骤6.2至步骤6.10，继续进行故障相关变量的选择，直至重构样本的T2统计量和SPE统计量均转为正常状态；
步骤7：重复上述步骤2至步骤6，直到f＝F，得到每种工业故障的故障相关变量集合；
步骤8：利用变量选择的结果，建立监测子空间及每个监测子空间的监测模型，具体步骤如下：
步骤8.1：将正常工况矩阵X分割为F+1个子空间为[Y1,Y2,…,Yf,…,YF+1]；
其中，当f<F+1时，
Figure BDA00019940548300000514
aj为矩阵X的第j个列向量也即为第j个可测得变量的正常工况样本集；当f＝F+1时，YF+1＝(aj)j∈{1,2,…,m}＝X；
步骤8.2：对于正常工况矩阵X的第f个子空间Yf，
若
Figure BDA00019940548300000515
则对第f个子空间中每个故障相关变量aj j∈Γf的样本集aj进行核密度估计，计算出每个故障相关变量aj的置信区间为95％的上限
Figure BDA0001994054830000061
若
Figure BDA0001994054830000062
则建立第f个子空间的PCA监测模型为
Figure BDA0001994054830000063
其中，
Figure BDA0001994054830000064
为矩阵Yf的负载矩阵，lf为矩阵Yf的主元个数，
Figure BDA0001994054830000065
为矩阵Yf的列数，
Figure BDA0001994054830000066
为矩阵Yf的单位矩阵，Ef为矩阵Yf的残差矩阵，计算第f个子空间的T2统计量的置信度为α的上限
Figure BDA0001994054830000067
和第f个子空间的SPE统计量的置信度为α的上限SPEfα；
步骤9：对工业生产过程进行工业故障的实时诊断，具体步骤如下：
步骤9.1：实时对工业生产过程中m个可测得变量的数据进行采集，得到初始待诊断样本，对初始待诊断样本进行标准化处理，得到待诊断矩阵xnew＝(a1,new,a2,new,…,aj,new,…,am,new)；其中，aj,new为待诊断矩阵xnew中第j个可测得变量的值；
步骤9.2：将待诊断矩阵xnew分割为F+1个子空间为[Y1,new,Y2,new,…,Yf,new,…,YF+1,new]；
其中，当f<F+1时，
Figure BDA0001994054830000068
当f＝F+1时，YF+1,new＝xnew；
步骤9.3：对于待诊断矩阵xnew的第f个子空间Yf,new，
若
Figure BDA0001994054830000069
则对第f个子空间Yf,new中每个故障相关变量aj,j∈Γf的值aj,new,j∈Γf与上限
Figure BDA00019940548300000610
进行比较：若有一个故障相关变量的值超过上限
Figure BDA00019940548300000611
则检测出有故障发生并诊断出故障类型为第f种故障，反之则检测出没有故障发生；
若
Figure BDA00019940548300000612
则计算第f个子空间Yf,new的样本
Figure BDA00019940548300000613
的T2统计量和SPE统计量分别为
Figure BDA00019940548300000614
Figure BDA00019940548300000615
其中，
Figure BDA00019940548300000616
为样本yf,new的得分向量，
Figure BDA00019940548300000617
将
Figure BDA00019940548300000618
与上限
Figure BDA00019940548300000619
与上限SPEfα进行比较：若
Figure BDA00019940548300000620
或SPEf,new>SPEfα，则检测出有故障发生，若f<F+1，则诊断出故障类型为第f种故障；若
Figure BDA00019940548300000621
且SPEf,new≤SPEfα，则检测出没有故障发生；
步骤9.4：重复上述步骤9.3，直到f＝F+1，完成每个子空间中的工业故障诊断。
所述步骤2中，正常工况样本的T2统计量的置信度为α的上限
Figure BDA0001994054830000071
为
Figure BDA0001994054830000072
其中，Fl,n-l；α代表具有l和n-l个自由度且置信水平为α的F分布的临界值；
正常工况样本的SPE统计量的置信度为α的上限SPEα为
Figure BDA0001994054830000073
其中，Cα是为高斯分布的(1-α)％的置信极限，
Figure BDA0001994054830000074
λ是协方差矩阵
Figure BDA0001994054830000075
的特征值。
本发明的有益效果为：
(1)本发明基于最小绝对值收敛及选择算子回归算法LASSO，同时考虑了正常工况数据与故障工况数据以及从正常工况数据到故障工况数据的相对变化，建立每种工业故障的故障相关变量选择模型，并基于结合故障重构思想的修正的最小角回归算法LARS对模型进行求解，相比于传统的仅采用正常工况数据进行建模、未充分利用历史故障工况数据信息且未考虑从正常工况数据到故障工况数据的相对变化的故障相关变量选择方法(贡献图方法)，能够更加精准地选择出每种故障的相关变量以进行故障根源分析。
(2)本发明基于故障相关变量选择的结果，基于故障子空间的思想建立PCA监测模型，实现工业故障的在线实时检测与诊断，相比于现有技术，能够在检测出工业故障的同时对故障类型进行诊断，提高工业故障诊断的效率及准确性。
附图说明
图1为青霉素发酵过程的流程图；
图2为本发明的基于相关变量选择的工业故障诊断方法的流程图；
图3为本发明的实施例一中正常工况样本的统计量示意图；
图4为本发明的实施例一中第1种工业故障的故障工况样本的统计量示意图；
图5为本发明的实施例一中第2种工业故障的故障工况样本的统计量示意图；
图6为统计量与故障相关变量及不相关变量的变化趋势对比图；
图7为本发明的实施例一中第1种工业故障的重构样本的统计量示意图；
图8为本发明的实施例一中第2种工业故障的重构样本的统计量示意图；
图9为本发明的实施例一中对工业故障的在线监测结果示意图；
图10为本发明的实施例二中对工业故障的在线监测结果示意图。
具体实施方式
下面将结合附图和具体实施方式，对本发明作进一步描述。
本发明的目的是提供一种基于相关变量选择的工业故障诊断方法，能够精准地提取每种工业故障的相关变量以进行故障根源分析，能够在检测出工业故障的同时对故障类型进行诊断，提高工业故障诊断的效率及准确性。
下面以青霉素发酵过程为例，对本发明作进一步描述。青霉素是世界各国需求量最大的抗生素，具有广泛的临床使用价值。在我国，经过多年的实践，关于青霉素发酵过程，业已积累了许多工业生产经验，但由于其化学结构复杂，初级代谢与次级代谢相互交叉，因此依据采集到的历史数据及相关领域专家的宝贵经验，进一步研究青霉素发酵过程的实时监测具有十分重要的意义。如图1所示，为青霉素发酵过程的流程图。整个青霉素生产过程主要包括4个生产时期：反映滞后期，菌体迅速生长期，青霉素合成期和菌体死亡期。首先，设置环境初始条件，进行微生物初始培养，然后再进行间歇补料发酵，不断生成青霉素。
实施例一
如图2所示，为本发明的基于相关变量选择的工业故障诊断方法的流程图。
本发明的基于相关变量选择的工业故障诊断方法，其特征在于，包括下述步骤：
步骤1：在工业生产过程中，确定m个可测得变量{a1,a2,…,am}，对正常工况下n个时刻的m个可测得变量的数据进行采集，得到n个正常工况样本，构成初始正常工况矩阵，对初始正常工况矩阵进行标准化处理，得到正常工况矩阵X＝(aij)n×m，aij为第i个正常工况样本中第j个可测得变量的值，i∈{1,2,…,n}，j∈{1,2,…,m}；确定F种工业故障，对每种工业故障下n个时刻的m个可测得变量的数据进行采集，得到每种工业故障的n个故障工况样本，构成每种工业故障的初始故障工况矩阵，对每种工业故障的初始故障工况矩阵进行标准化处理，得到第f种工业故障的故障工况矩阵Xf＝(aij f)n×m，f∈{1,2,...,F}，aij f为第f种工业故障的第i个故障工况样本中第j个可测得变量的值。
本实施例一中，首先确定青霉素生产过程的m＝16个可测得变量，如表1所示；确定F＝2种工业故障，如表2所示。在青霉素生产过程中，以1小时1次的采样频率，对正常工况下n＝400个时刻的16个可测得变量的数据进行采集，对每种工业故障下400个时刻的16个可测得变量的数据进行采集，得到400个正常工况样本及每种工业故障的400个故障工况样本，其中的部分样本如表3所示；通过标准化处理，得到正常工况矩阵X∈R400×16和故障工况矩阵X1∈R400×16和X2∈R400×16。
表1
Figure BDA0001994054830000091
表2
Figure BDA0001994054830000092
表3
Figure BDA0001994054830000093
步骤2：利用正常工况矩阵X建立PCA监测模型为
Figure BDA0001994054830000101
其中，P∈Rm×l为负载矩阵，l为矩阵X的主元个数，I∈Rm×m为单位矩阵，
Figure BDA0001994054830000102
为正常工况矩阵X在主元空间中的投影，E为正常工况矩阵X的残差矩阵；
计算每个正常工况样本x＝(ai1,ai2,…,aim)T的T2统计量和SPE统计量分别为
Figure BDA0001994054830000103
SPE＝||(I-PPT)x||2
其中，t＝PTx∈Rl×1为得分向量，
Figure BDA0001994054830000104
将n个正常工况样本的T2统计量和SPE统计量分别存入向量T2∈Rn×1和向量SPE∈Rn ×1中，并计算正常工况样本的T2统计量的置信度为α的上限
Figure BDA0001994054830000105
和SPE统计量的置信度为α的上限SPEα；其中，α＝95％。
在本实施例一中，所述步骤2中，正常工况样本的T2统计量的置信度为α的上限
Figure BDA0001994054830000106
为
Figure BDA0001994054830000107
其中，Fl,n-l；α代表具有l和n-l个自由度且置信水平为α的F分布的临界值；
正常工况样本的SPE统计量的置信度为α的上限SPEα为
Figure BDA0001994054830000108
其中，Cα是为高斯分布的(1-α)％的置信极限，
Figure BDA0001994054830000109
λ是协方差矩阵
Figure BDA00019940548300001010
的特征值。
步骤3：通过负载矩阵P对第f种工业故障的故障工况矩阵Xf进行分解：
Xf＝XfPPT+Xf(I-PPT)
计算第f种工业故障的每个故障工况样本的T2统计量和SPE统计量，并将第f种工业故障的n个故障工况样本的T2统计量和SPE统计量分别存入向量
Figure BDA00019940548300001011
和向量SPEf∈Rn×1中。
如图3所示，为本实施例一中正常工况样本的统计量示意图。从图3可以看出，正常工况样本的T2统计量基本位于上限
Figure BDA0001994054830000111
以下、SPE统计量基本位于上限SPEα以下。如图4、图5所示，为本发明的实施例一中第1种、第2种工业故障的故障工况样本的统计量示意图。从图4和图5可以看出，故障工况样本的T2统计量明显高于上限
Figure BDA0001994054830000112
SPE统计量明显高于上限SPEα。
步骤4：计算第f种工业故障的故障工况矩阵相对于正常工况矩阵的变化量矩阵为
ΔX＝Xf-X
计算第f种工业故障的向量
Figure BDA0001994054830000113
相对于向量T2的变化量矩阵为
Figure BDA0001994054830000114
计算第f种工业故障的向量SPEf相对于向量SPE的变化量矩阵为
ΔSPE＝SPEf-SPE
如图6所示，为统计量与故障相关变量及不相关变量的变化趋势对比图。从图6中可以看出，统计量的变化
Figure BDA00019940548300001113
和ΔSPE主要是由于ΔX中故障相关变量的变化引起的，从而只需比较统计量的波动与每个变量的变化趋势，找到其中最相关的若干变量即为故障相关变量。LASSO(最小绝对值收敛及选择算子，Least Absolute Shrinkage and SelectionOperator)利用L-1范数产生的稀疏效果来进行相应的变量选择，具体如下：
步骤5：基于最小绝对值收敛及选择算子回归算法LASSO，构建第f种工业故障的故障相关变量选择模型为
Figure BDA0001994054830000115
Figure BDA0001994054830000116
其中，
Figure BDA0001994054830000117
和βSPE分别为T2统计量和SPE统计量的回归系数向量，
Figure BDA0001994054830000118
和μSPE分别为T2统计量和SPE统计量的惩罚系数，
Figure BDA0001994054830000119
和μSPE分别用来控制
Figure BDA00019940548300001110
和βSPE的稀疏性。
由于LASSO问题不存在解析解，因此采用修正的最小角回归算法LARS(ModifiedLeast Angle Regression)来求解上述故障相关变量选择问题。LARS融合了故障重构的相关思想，可精确地确定故障相关变量的个数，具体如下：
步骤6：基于修正的最小角回归算法LARS，求解第f种工业故障的故障相关变量选择模型，具体步骤如下：
步骤6.1：对于第f种工业故障的故障相关变量选择模型的第一个子模型
Figure BDA00019940548300001111
确定第f种工业故障的在T2统计量下的估计向量
Figure BDA00019940548300001112
和第f种工业故障的在T2统计量下的故障相关变量集合
Figure BDA0001994054830000121
的初始值分别为
Figure BDA0001994054830000122
Figure BDA0001994054830000123
步骤6.2：计算相关性向量为
Figure BDA0001994054830000124
计算相关性向量
Figure BDA0001994054830000125
中元素的最大值为
Figure BDA0001994054830000126
将i添加至故障相关变量集合
Figure BDA0001994054830000127
中：
Figure BDA0001994054830000128
其中，
Figure BDA0001994054830000129
为相关性向量
Figure BDA00019940548300001210
的第i个元素，i∈{1,2,…,m}；
步骤6.3：令
Figure BDA00019940548300001211
构建矩阵
Figure BDA00019940548300001212
为
Figure BDA00019940548300001213
其中，
Figure BDA00019940548300001214
为
Figure BDA00019940548300001215
的符号，
Figure BDA00019940548300001216
为矩阵ΔX的第i列向量，
Figure BDA00019940548300001217
为全1向量，
Figure BDA00019940548300001218
为故障相关变量集合
Figure BDA00019940548300001219
中的元素个数；
步骤6.4：计算估计向量
Figure BDA00019940548300001220
的更新方向
Figure BDA00019940548300001221
和更新步长
Figure BDA00019940548300001222
分别为：
Figure BDA00019940548300001223
Figure BDA00019940548300001224
其中，
Figure BDA00019940548300001225
为故障相关变量集合
Figure BDA00019940548300001226
的补集，aj为向量
Figure BDA00019940548300001227
中的第j个元素；
步骤6.5：更新估计向量
Figure BDA00019940548300001228
步骤6.6：采用上述步骤6.1至步骤6.5的方法，对第f种工业故障的故障相关变量选择模型的第二个子模型
Figure BDA00019940548300001229
进行求解，得到第f种工业故障的在SPE统计量下的估计向量
Figure BDA00019940548300001230
和第f种工业故障的在SPE统计量下的故障相关变量集合ΓSPE；
步骤6.7：对集合
Figure BDA00019940548300001231
与ΓSPE求取并集，得到第f种工业故障的故障相关变量集合为
Figure BDA00019940548300001232
步骤6.8：利用故障相关变量集合Γf和T2统计量，对第f种工业故障的每个故障工况样本xf＝(ai1 f,ai2 f,…,aim f)T进行重构，建立第一重构模型为
Figure BDA0001994054830000131
其中，
Figure BDA0001994054830000132
为故障方向矩阵，
Figure BDA0001994054830000133
为第f种工业故障的故障相关变量集合Γf中元素的个数也即故障相关变量的个数，矩阵Ξ的第i行第j列的元素为
Figure BDA0001994054830000134
(Γf)j为集合Γf中的第j个元素；
Figure BDA0001994054830000135
为第一幅值估计向量；
求解第一重构模型，得到当前最优的第一幅值估计向量为
Figure BDA0001994054830000136
从而，得到第f种工业故障的第一重构样本为
Figure BDA0001994054830000137
步骤6.9：利用故障相关变量集合Γf和SPE统计量，对第f种工业故障的每个故障工况样本xf＝(ai1 f,ai2 f,…,aim f)T进行重构，建立第二重构模型为
Figure BDA0001994054830000138
其中，
Figure BDA0001994054830000139
为第二幅值估计向量；
求解第二重构模型，得到当前最优的第二幅值估计向量为
hSPE＝(ΞT(I-PPT)Ξ)-1ΞT(I-PPT)xf
从而，得到第f种工业故障的第二重构样本为
zSPE＝xf-ΞhSPE
＝xf-Ξ(ΞT(I-PPT)Ξ)-1ΞT(I-PPT)xf
步骤6.10：计算第f种工业故障的每个第一重构样本
Figure BDA00019940548300001310
的T2统计量为
Figure BDA00019940548300001311
其中，
Figure BDA00019940548300001312
为样本
Figure BDA00019940548300001313
的得分向量；
计算第f种工业故障的每个第二重构样本zSPE∈Rm×1的SPE统计量为
Figure BDA0001994054830000141
将第f种工业故障的n个第一重构样本的T2统计量和n个第二重构样本的SPE统计量分别存入向量
Figure BDA0001994054830000142
和向量
Figure BDA0001994054830000143
中；
步骤6.11：将向量
Figure BDA0001994054830000144
中的每个元素与上限
Figure BDA0001994054830000145
进行比较、将向量
Figure BDA0001994054830000146
中的每个元素与上限SPEα进行比较：
若向量
Figure BDA0001994054830000147
中有98％以上的元素低于上限
Figure BDA0001994054830000148
且向量
Figure BDA0001994054830000149
中有98％以上的元素低于上限SPEα，则认为重构样本的T2统计量和SPE统计量均转为正常状态，从而第f种工业故障的所有故障相关变量都已经位于集合Γf中，进入步骤7；
否则，重复上述步骤6.2至步骤6.10，继续进行故障相关变量的选择，直至重构样本的T2统计量和SPE统计量均转为正常状态。
步骤7：重复上述步骤2至步骤6，直到f＝F，得到每种工业故障的故障相关变量集合。
本实施例一中，当f＝1也即对于第一种工业故障，在一次迭代过程中，步骤6.7中得到第1种工业故障的故障相关变量集合为Γ1＝{1}，在步骤6.8、步骤6.9中对第1种工业故障的每个故障工况样本进行重构，得到第1种工业故障的第一重构样本和第二重构样本，在步骤6.10中计算得到第1种工业故障的每个第一重构样本的T2统计量和每个第二重构样本的SPE统计量；在步骤6.11中将向量
Figure BDA00019940548300001410
中的每个元素与上限
Figure BDA00019940548300001411
进行比较、将向量
Figure BDA00019940548300001412
中的每个元素与上限SPEα进行比较，发现向量
Figure BDA00019940548300001413
中有98％以上的元素低于上限
Figure BDA00019940548300001414
且向量
Figure BDA00019940548300001415
中有98％以上的元素低于上限SPEα，从而重构样本的T2统计量和SPE统计量均转为正常状态，第1种工业故障的所有故障相关变量都已经位于集合Γf中，从而确定第1种工业故障的故障相关变量为a1。重复步骤2至步骤6，对第2种工业故障的故障相关变量进行选择，最终得到第2种工业故障的故障相关变量集合为Γ2＝{2}，从而第2种工业故障的故障相关变量为a2。
如图7、图8所示，分别为本发明的实施例一中第1种、第2种工业故障的重构样本的统计量示意图。从图7和图8可以看出，经过重构后的样本的统计量与正常样本的统计量非常相似，证明故障变量影响已被消除，过程恢复到正常状态。
步骤8：利用变量选择的结果，建立监测子空间及每个监测子空间的监测模型，具体步骤如下：
步骤8.1：将正常工况矩阵X分割为F+1个子空间为[Y1,Y2,…,Yf,…,YF+1]；
其中，当f<F+1时，
Figure BDA00019940548300001416
aj为矩阵X的第j个列向量也即为第j个可测得变量的正常工况样本集；当f＝F+1时，YF+1＝(aj)j∈{1,2,…,m}＝X；
步骤8.2：对于正常工况矩阵X的第f个子空间Yf，
若
Figure BDA0001994054830000151
则对第f个子空间中每个故障相关变量aj j∈Γf的样本集aj进行核密度估计，计算出每个故障相关变量aj的置信区间为95％的上限
Figure BDA0001994054830000152
若
Figure BDA0001994054830000153
则建立第f个子空间的PCA监测模型为
Figure BDA0001994054830000154
其中，
Figure BDA0001994054830000155
为矩阵Yf的负载矩阵，lf为矩阵Yf的主元个数，
Figure BDA0001994054830000156
为矩阵Yf的列数，
Figure BDA0001994054830000157
为矩阵Yf的单位矩阵，Ef为矩阵Yf的残差矩阵，计算第f个子空间的T2统计量的置信度为α的上限
Figure BDA0001994054830000158
和第f个子空间的SPE统计量的置信度为α的上限SPEfα；
本实施例一中，将正常工况矩阵X分割为3个子空间[Y1,Y2,Y3]，Y1＝(a1)∈Rn×1，Y2＝(a2)∈Rn×1，Y3＝X。对于第1个和第2个子空间，其包含的变量个数均小于3，也即
Figure BDA0001994054830000159
n2<3，从而对第1个子空间中故障相关变量a1的样本集a1进行核密度估计，计算出第1个子空间中故障相关变量a1的置信区间为95％的上限
Figure BDA00019940548300001510
对第2个子空间中故障相关变量a2的样本集a2进行核密度估计，计算出第2个子空间中故障相关变量a2的置信区间为95％的上限
Figure BDA00019940548300001511
对于第3个子空间，其包含的变量个数大于3，也即
Figure BDA00019940548300001512
从而建立第3个子空间的PCA监测模型，并计算第3个子空间的T2统计量的置信度为α的上限
Figure BDA00019940548300001513
和第3个子空间的SPE统计量的置信度为α的上限SPE3α。
步骤9：对工业生产过程进行工业故障的实时诊断，具体步骤如下：
步骤9.1：实时对工业生产过程中m个可测得变量的数据进行采集，得到初始待诊断样本，对初始待诊断样本进行标准化处理，得到待诊断矩阵xnew＝(a1,new,a2,new,…,aj,new,…,am,new)；其中，aj,new为待诊断矩阵xnew中第j个可测得变量的值；
步骤9.2：将待诊断矩阵xnew分割为F+1个子空间为[Y1,new,Y2,new,…,Yf,new,…,YF+1,new]；
其中，当f<F+1时，
Figure BDA00019940548300001514
当f＝F+1时，YF+1,new＝xnew；
步骤9.3：对于待诊断矩阵xnew的第f个子空间Yf,new，
若
Figure BDA0001994054830000161
则对第f个子空间Yf,new中每个故障相关变量aj,j∈Γf的值aj,new,j∈Γf与上限
Figure BDA0001994054830000162
进行比较：若有一个故障相关变量的值超过上限
Figure BDA0001994054830000163
则检测出有故障发生并诊断出故障类型为第f种故障，反之则检测出没有故障发生；
若
Figure BDA0001994054830000164
则计算第f个子空间Yf,new的样本
Figure BDA0001994054830000165
的T2统计量和SPE统计量分别为
Figure BDA0001994054830000166
Figure BDA0001994054830000167
其中，
Figure BDA0001994054830000168
为样本yf,new的得分向量，
Figure BDA0001994054830000169
将
Figure BDA00019940548300001610
与上限
Figure BDA00019940548300001611
与上限SPEfα进行比较：若
Figure BDA00019940548300001612
或SPEf,new>SPEfα，则检测出有故障发生，若f<F+1，则诊断出故障类型为第f种故障；若
Figure BDA00019940548300001613
且SPEf,new≤SPEfα，则检测出没有故障发生；
步骤9.4：重复上述步骤9.3，直到f＝F+1，完成每个子空间中的工业故障诊断。
本实施例一中，实时对青霉素生产过程中16个可测得变量的数据进行采集，经过标准化处理后得到待诊断矩阵xnew，将矩阵xnew分割为3个子空间为[Y1,new,Y2,new,Y3,new]，Y1,new＝(a1,new)，Y2,new＝(a2,new)，Y3,new＝xnew。在第1个子空间Y1,new中，比较故障相关变量a1的值a1,new与上限
Figure BDA00019940548300001614
在第2个子空间Y2,new中，比较故障相关变量a2的值a2,new与上限
Figure BDA00019940548300001615
在第3个子空间Y3,new中，计算样本
Figure BDA00019940548300001616
的T2统计量和SPE统计量分别为
Figure BDA00019940548300001617
并比较
Figure BDA00019940548300001618
与上限
Figure BDA00019940548300001619
与上限SPE3α。
如图9所示，为本发明的实施例一中对工业故障的在线监测结果示意图。在本实施例一中，实时采集400个时刻的待诊断样本。从图9可以看出，在子空间1中，在第201个样本处，变量a1的值超过上限
Figure BDA00019940548300001620
在子空间2中，变量a2的值始终在上限
Figure BDA00019940548300001621
之下；在子空间3中，在第201个样本处，样本y3,new的T2统计量和SPE统计量均分别超过上限
Figure BDA00019940548300001622
和上限SPE3α。从而，在子空间1和子空间3中均检测出有故障发生，并在子空间1中诊断出故障类型为第1种故障，在子空间2中则检测出没有故障发生。从而，本实施例一中，通过本发明的方法，检测出在第201个时刻有故障发生并诊断出故障类型为第1种故障。
实施例二
如图10所示，为本发明的实施例二中对工业故障的在线监测结果示意图。本实施例二与上述实施例一的不同之处在于：在子空间1中，变量a1的值始终在上限
Figure BDA0001994054830000171
之下；在子空间2中，在第201个样本处，变量a2的值超过上限
Figure BDA0001994054830000172
在子空间3中，在第201个样本处，样本y3,new的T2统计量和SPE统计量均分别超过上限
Figure BDA0001994054830000173
和上限SPE3α。从而，在子空间2和子空间3中均检测出有故障发生，并在子空间2中诊断出故障类型为第2种故障，在子空间1中则检测出没有故障发生。从而，本实施例二中，通过本发明的方法，检测出在第201个时刻有故障发生并诊断出故障类型为第2种故障。
显然，上述实施例仅仅是本发明的一部分实施例，而不是全部的实施例。上述实施例仅用于解释本发明，并不构成对本发明保护范围的限定。基于上述实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，也即凡在本申请的精神和原理之内所作的所有修改、等同替换和改进等，均落在本发明要求的保护范围内。

Claims (2)

1.一种基于相关变量选择的工业故障诊断方法，其特征在于，包括下述步骤：
步骤1：在工业生产过程中，确定m个可测得变量{a1,a2,…,am}，对正常工况下n个时刻的m个可测得变量的数据进行采集，得到n个正常工况样本，构成初始正常工况矩阵，对初始正常工况矩阵进行标准化处理，得到正常工况矩阵X＝(aij)n×m，aij为第i个正常工况样本中第j个可测得变量的值，i∈{1,2,…,n}，j∈{1,2,…,m}；确定F种工业故障，对每种工业故障下n个时刻的m个可测得变量的数据进行采集，得到每种工业故障的n个故障工况样本，构成每种工业故障的初始故障工况矩阵，对每种工业故障的初始故障工况矩阵进行标准化处理，得到第f种工业故障的故障工况矩阵Xf＝(aij f)n×m，f∈{1,2,...,F}，aij f为第f种工业故障的第i个故障工况样本中第j个可测得变量的值；
步骤2：利用正常工况矩阵X建立PCA监测模型为
Figure FDA0001994054820000011
其中，P∈Rm×l为负载矩阵，l为矩阵X的主元个数，I∈Rm×m为单位矩阵，
Figure FDA0001994054820000012
为正常工况矩阵X在主元空间中的投影，E为正常工况矩阵X的残差矩阵；
计算每个正常工况样本x＝(ai1,ai2,…,aim)T的T2统计量和SPE统计量分别为
Figure FDA0001994054820000013
SPE＝||(I-PPT)x||2
其中，t＝PTx∈Rl×1为得分向量，
Figure FDA0001994054820000014
将n个正常工况样本的T2统计量和SPE统计量分别存入向量T2∈Rn×1和向量SPE∈Rn×1中，并计算正常工况样本的T2统计量的置信度为α的上限
Figure FDA0001994054820000015
和SPE统计量的置信度为α的上限SPEα；其中，α＝95％；
步骤3：通过负载矩阵P对第f种工业故障的故障工况矩阵Xf进行分解：
Xf＝XfPPT+Xf(I-PPT)
计算第f种工业故障的每个故障工况样本的T2统计量和SPE统计量，并将第f种工业故障的n个故障工况样本的T2统计量和SPE统计量分别存入向量
Figure FDA0001994054820000016
和向量SPEf∈Rn×1中；
步骤4：计算第f种工业故障的故障工况矩阵相对于正常工况矩阵的变化量矩阵为
ΔX＝Xf-X
计算第f种工业故障的向量
Figure FDA0001994054820000021
相对于向量T2的变化量矩阵为
Figure FDA0001994054820000022
计算第f种工业故障的向量SPEf相对于向量SPE的变化量矩阵为
ΔSPE＝SPEf-SPE
步骤5：基于最小绝对值收敛及选择算子回归算法LASSO，构建第f种工业故障的故障相关变量选择模型为
Figure FDA0001994054820000023
Figure FDA0001994054820000024
其中，βT2和βSPE分别为T2统计量和SPE统计量的回归系数向量，
Figure FDA0001994054820000025
和μSPE分别为T2统计量和SPE统计量的惩罚系数，
Figure FDA0001994054820000026
和μSPE分别用来控制
Figure FDA0001994054820000027
和βSPE的稀疏性；
步骤6：基于修正的最小角回归算法LARS，求解第f种工业故障的故障相关变量选择模型，具体步骤如下：
步骤6.1：对于第f种工业故障的故障相关变量选择模型的第一个子模型
Figure FDA0001994054820000028
确定第f种工业故障的在T2统计量下的估计向量
Figure FDA0001994054820000029
和第f种工业故障的在T2统计量下的故障相关变量集合
Figure FDA00019940548200000210
的初始值分别为
Figure FDA00019940548200000211
Figure FDA00019940548200000212
步骤6.2：计算相关性向量为
Figure FDA00019940548200000213
计算相关性向量
Figure FDA00019940548200000214
中元素的最大值为
Figure FDA00019940548200000215
将i添加至故障相关变量集合
Figure FDA00019940548200000216
中：
Figure FDA00019940548200000217
其中，
Figure FDA00019940548200000218
为相关性向量
Figure FDA00019940548200000219
的第i个元素，i∈{1,2,…,m}；
步骤6.3：令
Figure FDA00019940548200000220
构建矩阵
Figure FDA0001994054820000031
为
Figure FDA0001994054820000032
其中，
Figure FDA0001994054820000033
为
Figure FDA0001994054820000034
的符号，
Figure FDA0001994054820000035
为矩阵ΔX的第i列向量，
Figure FDA0001994054820000036
为全1向量，
Figure FDA0001994054820000037
为故障相关变量集合
Figure FDA0001994054820000038
中的元素个数；
步骤6.4：计算估计向量
Figure FDA0001994054820000039
的更新方向
Figure FDA00019940548200000310
和更新步长
Figure FDA00019940548200000311
分别为：
Figure FDA00019940548200000312
Figure FDA00019940548200000313
其中，
Figure FDA00019940548200000314
为故障相关变量集合
Figure FDA00019940548200000315
的补集，aj为向量
Figure FDA00019940548200000316
中的第j个元素；
步骤6.5：更新估计向量
Figure FDA00019940548200000317
步骤6.6：采用上述步骤6.1至步骤6.5的方法，对第f种工业故障的故障相关变量选择模型的第二个子模型
Figure FDA00019940548200000318
进行求解，得到第f种工业故障的在SPE统计量下的估计向量
Figure FDA00019940548200000319
和第f种工业故障的在SPE统计量下的故障相关变量集合ΓSPE；
步骤6.7：对集合
Figure FDA00019940548200000320
与ΓSPE求取并集，得到第f种工业故障的故障相关变量集合为
Figure FDA00019940548200000321
步骤6.8：利用故障相关变量集合Γf和T2统计量，对第f种工业故障的每个故障工况样本xf＝(ai1 f,ai2 f,…,aim f)T进行重构，建立第一重构模型为
Figure FDA00019940548200000322
其中，
Figure FDA00019940548200000323
为故障方向矩阵，
Figure FDA00019940548200000324
为第f种工业故障的故障相关变量集合Γf中元素的个数也即故障相关变量的个数，矩阵Ξ的第i行第j列的元素为
Figure FDA00019940548200000325
(Γf)j为集合Γf中的第j个元素；
Figure FDA00019940548200000326
为第一幅值估计向量；
求解第一重构模型，得到当前最优的第一幅值估计向量为
Figure FDA00019940548200000327
从而，得到第f种工业故障的第一重构样本为
Figure FDA0001994054820000041
步骤6.9：利用故障相关变量集合Γf和SPE统计量，对第f种工业故障的每个故障工况样本xf＝(ai1 f,ai2 f,…,aim f)T进行重构，建立第二重构模型为
Figure FDA0001994054820000042
其中，
Figure FDA0001994054820000043
为第二幅值估计向量；
求解第二重构模型，得到当前最优的第二幅值估计向量为
hSPE＝(ΞT(I-PPT)Ξ)-1ΞT(I-PPT)xf
从而，得到第f种工业故障的第二重构样本为
zSPE＝xf-ΞhSPE
＝xf-Ξ(ΞT(I-PPT)Ξ)-1ΞT(I-PPT)xf
步骤6.10：计算第f种工业故障的每个第一重构样本
Figure FDA0001994054820000044
的T2统计量为
Figure FDA0001994054820000045
其中，
Figure FDA0001994054820000046
为样本
Figure FDA0001994054820000047
的得分向量；
计算第f种工业故障的每个第二重构样本zSPE∈Rm×1的SPE统计量为
Figure FDA0001994054820000048
将第f种工业故障的n个第一重构样本的T2统计量和n个第二重构样本的SPE统计量分别存入向量
Figure FDA0001994054820000049
和向量
Figure FDA00019940548200000410
中；
步骤6.11：将向量
Figure FDA00019940548200000411
中的每个元素与上限
Figure FDA00019940548200000412
进行比较、将向量
Figure FDA00019940548200000413
中的每个元素与上限SPEα进行比较：
若向量
Figure FDA00019940548200000414
中有98％以上的元素低于上限
Figure FDA00019940548200000415
且向量
Figure FDA00019940548200000416
中有98％以上的元素低于上限SPEα，则认为重构样本的T2统计量和SPE统计量均转为正常状态，从而第f种工业故障的所有故障相关变量都已经位于集合Γf中，进入步骤7；
否则，重复上述步骤6.2至步骤6.10，继续进行故障相关变量的选择，直至重构样本的T2统计量和SPE统计量均转为正常状态；
步骤7：重复上述步骤2至步骤6，直到f＝F，得到每种工业故障的故障相关变量集合；
步骤8：利用变量选择的结果，建立监测子空间及每个监测子空间的监测模型，具体步骤如下：
步骤8.1：将正常工况矩阵X分割为F+1个子空间为[Y1,Y2,…,Yf,…,YF+1]；
其中，当f<F+1时，
Figure FDA0001994054820000051
aj为矩阵X的第j个列向量也即为第j个可测得变量的正常工况样本集；当f＝F+1时，YF+1＝(aj)j∈{1,2,…,m}＝X；
步骤8.2：对于正常工况矩阵X的第f个子空间Yf，
若
Figure FDA0001994054820000052
则对第f个子空间中每个故障相关变量aj j∈Γf的样本集aj进行核密度估计，计算出每个故障相关变量aj的置信区间为95％的上限
Figure FDA0001994054820000053
若
Figure FDA0001994054820000054
则建立第f个子空间的PCA监测模型为
Figure FDA0001994054820000055
其中，
Figure FDA0001994054820000056
为矩阵Yf的负载矩阵，lf为矩阵Yf的主元个数，
Figure FDA0001994054820000057
为矩阵Yf的列数，
Figure FDA0001994054820000058
为矩阵Yf的单位矩阵，Ef为矩阵Yf的残差矩阵，计算第f个子空间的T2统计量的置信度为α的上限
Figure FDA0001994054820000059
和第f个子空间的SPE统计量的置信度为α的上限SPEfα；
步骤9：对工业生产过程进行工业故障的实时诊断，具体步骤如下：
步骤9.1：实时对工业生产过程中m个可测得变量的数据进行采集，得到初始待诊断样本，对初始待诊断样本进行标准化处理，得到待诊断矩阵xnew＝(a1,new,a2,new,…,aj,new,…,am,new)；其中，aj,new为待诊断矩阵xnew中第j个可测得变量的值；
步骤9.2：将待诊断矩阵xnew分割为F+1个子空间为[Y1,new,Y2,new,…,Yf,new,…,YF+1,new]；
其中，当f<F+1时，
Figure FDA00019940548200000510
当f＝F+1时，YF+1,new＝xnew；
步骤9.3：对于待诊断矩阵xnew的第f个子空间Yf,new，
若
Figure FDA00019940548200000511
则对第f个子空间Yf,new中每个故障相关变量aj,j∈Γf的值aj,new,j∈Γf与上限
Figure FDA00019940548200000512
进行比较：若有一个故障相关变量的值超过上限
Figure FDA00019940548200000513
则检测出有故障发生并诊断出故障类型为第f种故障，反之则检测出没有故障发生；
若
Figure FDA00019940548200000514
则计算第f个子空间Yf,new的样本
Figure FDA00019940548200000515
的T2统计量和SPE统计量分别为
Figure FDA0001994054820000061
Figure FDA0001994054820000062
其中，
Figure FDA0001994054820000063
为样本yf,new的得分向量，
Figure FDA0001994054820000064
将
Figure FDA0001994054820000065
与上限
Figure FDA0001994054820000066
SPEf,new与上限SPEfα进行比较：若
Figure FDA0001994054820000067
或SPEf,new>SPEfα，则检测出有故障发生，若f<F+1，则诊断出故障类型为第f种故障；若
Figure FDA0001994054820000068
且SPEf,new≤SPEfα，则检测出没有故障发生；
步骤9.4：重复上述步骤9.3，直到f＝F+1，完成每个子空间中的工业故障诊断。
2.根据权利要求1所述的基于相关变量选择的工业故障诊断方法，其特征在于，所述步骤2中，正常工况样本的T2统计量的置信度为α的上限
Figure FDA0001994054820000069
为
Figure FDA00019940548200000610
其中，Fl,n-l；α代表具有l和n-l个自由度且置信水平为α的F分布的临界值；
正常工况样本的SPE统计量的置信度为α的上限SPEα为
Figure FDA00019940548200000611
其中，Cα是为高斯分布的(1-α)％的置信极限，
Figure FDA00019940548200000612
λ是协方差矩阵
Figure FDA00019940548200000613
的特征值。
CN201910189841.4A 2019-03-13 2019-03-13 一种基于相关变量选择的工业故障诊断方法 Active CN110083860B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910189841.4A CN110083860B (zh) 2019-03-13 2019-03-13 一种基于相关变量选择的工业故障诊断方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910189841.4A CN110083860B (zh) 2019-03-13 2019-03-13 一种基于相关变量选择的工业故障诊断方法

Publications (2)

Publication Number Publication Date
CN110083860A CN110083860A (zh) 2019-08-02
CN110083860B true CN110083860B (zh) 2023-01-13

Family

ID=67413209

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910189841.4A Active CN110083860B (zh) 2019-03-13 2019-03-13 一种基于相关变量选择的工业故障诊断方法

Country Status (1)

Country Link
CN (1) CN110083860B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111651220B (zh) * 2020-06-04 2023-08-18 上海电力大学 一种基于深度强化学习的Spark参数自动优化方法及系统
CN113076211B (zh) * 2021-03-29 2024-02-23 中国人民解放军火箭军工程大学 一种基于故障重构的质量相关故障诊断及误报警反馈方法
CN113189968B (zh) * 2021-05-08 2022-08-26 哈尔滨工业大学 互联工业过程的分布式故障诊断方法
CN113359679A (zh) * 2021-06-24 2021-09-07 东北大学 基于重构幅值趋势特征的工业过程故障诊断方法
CN113554061B (zh) * 2021-06-25 2022-11-22 东南大学 重构pca算法中主元个数的选择方法
CN116300774B (zh) * 2023-05-23 2023-08-08 蓝星智云(山东)智能科技有限公司 基于主元分析与核密度估计的间歇过程可视化监控方法
CN116679669B (zh) * 2023-06-07 2024-03-26 矿冶科技集团有限公司 一种筛分系统故障诊断方法及系统

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
AR092747A1 (es) * 2013-09-30 2015-04-29 Ypf Tecnologia Sa Dispositivo y procedimiento para deteccion y/o diagnostico de fallas en procesos, equipos y sensores
CN103729562A (zh) * 2013-12-31 2014-04-16 东北大学 一种基于重构识别分析的青霉素发酵过程故障监测方法
CN109062189B (zh) * 2018-08-30 2020-06-30 华中科技大学 一种用于复杂故障的工业过程故障诊断方法

Also Published As

Publication number Publication date
CN110083860A (zh) 2019-08-02

Similar Documents

Publication Publication Date Title
CN110083860B (zh) 一种基于相关变量选择的工业故障诊断方法
CN110336534B (zh) 一种基于光伏阵列电气参数时间序列特征提取的故障诊断方法
CN108803520B (zh) 一种基于变量非线性自相关性剔除的动态过程监测方法
CN107632592B (zh) 基于高效递推核主元分析的非线性时变过程故障监测方法
CN108897286B (zh) 一种基于分散式非线性动态关系模型的故障检测方法
CN101458522A (zh) 基于主元分析和支持向量数据描述的多工况过程监控方法
CN111275288A (zh) 基于XGBoost的多维数据异常检测方法与装置
CN108549908B (zh) 基于多采样概率核主成分模型的化工过程故障检测方法
CN112817291B (zh) 基于混合特性评价和子空间分解的层次故障监测方法
CN111340110A (zh) 一种基于工业过程运行状态趋势分析的故障预警方法
CN110032799A (zh) 一种微生物制药过程的角相似度阶段划分及监测方法
CN116383636A (zh) 一种基于pca与lstm融合算法的磨煤机故障预警方法
CN110942258B (zh) 一种性能驱动的工业过程异常监测方法
CN110490496B (zh) 一种基于分步约简筛选复杂工业过程中影响产品质量的敏感变量的方法
CN113703422B (zh) 一种基于特征分析处理的燃气轮机气动执行机构故障诊断方法
CN116401516A (zh) 一种基于深度学习的电力负荷异常数据检测与修正方法
CN116431966A (zh) 一种增量式特征解耦自编码器的堆芯温度异常检测方法
CN114354184A (zh) 一种基于深度学习的大型回转装备主轴健康预警模型建立方法和装置
CN111188761B (zh) 一种基于Fourier-CVA模型面向机泵设备的监测方法
CN112149054B (zh) 基于时序扩展的正交邻域保持嵌入模型的构建与应用
CN111914886B (zh) 一种基于在线简略核学习的非线性化工过程监测方法
CN111914471A (zh) 一种基于快速核独立成分分析的精馏塔故障检测方法
CN109522657B (zh) 一种基于相关性网络和svdd的燃气轮机异常检测方法
CN116339275A (zh) 基于全结构动态自回归隐变量模型的多尺度过程故障检测方法
CN114707424A (zh) 基于质量相关慢特征分析算法的化工过程软测量方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant