CN109143995A - 一种基于质量相关慢特征充分分解的闭环***精细运行状态监测方法 - Google Patents
一种基于质量相关慢特征充分分解的闭环***精细运行状态监测方法 Download PDFInfo
- Publication number
- CN109143995A CN109143995A CN201810771664.6A CN201810771664A CN109143995A CN 109143995 A CN109143995 A CN 109143995A CN 201810771664 A CN201810771664 A CN 201810771664A CN 109143995 A CN109143995 A CN 109143995A
- Authority
- CN
- China
- Prior art keywords
- matrix
- quality
- variable
- space
- feature
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 100
- 238000012544 monitoring process Methods 0.000 title claims abstract description 32
- 230000008569 process Effects 0.000 claims abstract description 62
- 230000003068 static effect Effects 0.000 claims abstract description 62
- 230000002159 abnormal effect Effects 0.000 claims abstract description 11
- 239000011159 matrix material Substances 0.000 claims description 122
- 230000008859 change Effects 0.000 claims description 15
- 238000004458 analytical method Methods 0.000 claims description 13
- 238000000354 decomposition reaction Methods 0.000 claims description 11
- 238000004519 manufacturing process Methods 0.000 claims description 8
- 238000005259 measurement Methods 0.000 claims description 5
- 238000010219 correlation analysis Methods 0.000 claims description 4
- 238000005457 optimization Methods 0.000 claims description 4
- 230000009471 action Effects 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 238000004540 process dynamic Methods 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims 1
- 238000001514 detection method Methods 0.000 abstract description 8
- 238000000605 extraction Methods 0.000 abstract description 3
- 239000000047 product Substances 0.000 description 41
- 238000011160 research Methods 0.000 description 8
- 238000010586 diagram Methods 0.000 description 7
- 230000004888 barrier function Effects 0.000 description 4
- 238000003745 diagnosis Methods 0.000 description 3
- 238000012423 maintenance Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 2
- 230000004069 differentiation Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 230000005856 abnormality Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 239000006227 byproduct Substances 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000012824 chemical production Methods 0.000 description 1
- 238000011217 control strategy Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- MXCPYJZDGPQDRA-UHFFFAOYSA-N dialuminum;2-acetyloxybenzoic acid;oxygen(2-) Chemical compound [O-2].[O-2].[O-2].[Al+3].[Al+3].CC(=O)OC1=CC=CC=C1C(O)=O MXCPYJZDGPQDRA-UHFFFAOYSA-N 0.000 description 1
- 239000007789 gas Substances 0.000 description 1
- 239000011261 inert gas Substances 0.000 description 1
- 239000012263 liquid product Substances 0.000 description 1
- 230000007257 malfunction Effects 0.000 description 1
- 239000002994 raw material Substances 0.000 description 1
- 230000008439 repair process Effects 0.000 description 1
- 238000012827 research and development Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 230000026676 system process Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B19/00—Programme-control systems
- G05B19/02—Programme-control systems electric
- G05B19/418—Total factory control, i.e. centrally controlling a plurality of machines, e.g. direct or distributed numerical control [DNC], flexible manufacturing systems [FMS], integrated manufacturing systems [IMS] or computer integrated manufacturing [CIM]
- G05B19/41875—Total factory control, i.e. centrally controlling a plurality of machines, e.g. direct or distributed numerical control [DNC], flexible manufacturing systems [FMS], integrated manufacturing systems [IMS] or computer integrated manufacturing [CIM] characterised by quality surveillance of production
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
- G05B19/00—Programme-control systems
- G05B19/02—Programme-control systems electric
- G05B19/418—Total factory control, i.e. centrally controlling a plurality of machines, e.g. direct or distributed numerical control [DNC], flexible manufacturing systems [FMS], integrated manufacturing systems [IMS] or computer integrated manufacturing [CIM]
- G05B19/41885—Total factory control, i.e. centrally controlling a plurality of machines, e.g. direct or distributed numerical control [DNC], flexible manufacturing systems [FMS], integrated manufacturing systems [IMS] or computer integrated manufacturing [CIM] characterised by modeling, simulation of the manufacturing system
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
- Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
- Y02P90/02—Total factory control, e.g. smart factories, flexible manufacturing systems [FMS] or integrated manufacturing systems [IMS]
Landscapes
- Engineering & Computer Science (AREA)
- Manufacturing & Machinery (AREA)
- General Engineering & Computer Science (AREA)
- Quality & Reliability (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Automation & Control Theory (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Feedback Control In General (AREA)
Abstract
本发明公开了一种基于质量相关慢特征充分分解的闭环***精细运行状态监测方法。本方法通过提取过程变量与质量变量的静态相关性，将闭环大规模过程变量分解成与质量相关和与质量无关的两个空间，并在每个空间中提取过程的静态和动态信息对其进行建模，分别在两个空间内建立静态指标和动态指标，针对两个空间进行协同的故障检测。该方法不仅可以实现对闭环***动态异常情况和工况变化的监测，有效识别闭环故障***中真正故障的发生，并且通过对变量空间进行质量相关分解，可以进一步实现对闭环***产品质量的相关监测，判断出故障和异常具体发生在哪个空间，是否影响产品质量。
Description
技术领域
本发明属于工业***过程监测领域，特别是针对闭环***产品质量相关的精细运行状态监测方法。
背景技术
随着科学技术的发展，现代工业***的规模和复杂程度都在日益提高。复杂大***一旦出现异常，都可能带来重大的财产损失和人员伤亡。通过对***设备日常的运行状态监测，提前检测出设备异常，使设备得到及时的维修，不仅使得检修周期延长，维修费用大大降低，而且检修更有针对性，设备使用寿命得到延长，避免了严重安全事故的发生。因此，为了保证***在运行的安全性与可靠性，及时发现***运行中的异常情况并进行处理，减少生产中的安全隐患，提高设备使用周期，十分有必要采用有效手段对***进行实时监测和故障检测。同时随着传感技术的发展，在工业现场获得数据变得越来越容易，过程数据中蕴含了大量的过程信息，基于数据的状态监测和故障监测逐渐成为研究的热点。
在过去的几十年时间里，过程监测和故障检测技术得到了广泛研究和发展，大量研究成果得到发表，前人对基于数据的故障检测和故障诊断作出了相应的研究。主成分分析(PCA)，偏最小二乘(PLS)和费舍尔判别分析(FDA)等多元统计分析方法已经被广泛应用于基于数据的过程监测领域。然而现有的大部分的过程监测和故障检测方法都是针对开环***设计的，没有考虑到闭环控制规律的影响。但是在真实工业***中，为满足稳定性、快速性、准确性要求，达到实际生产目标保证产品质量维持在合理的水平，需要对***施加闭环反馈控制，大量闭环控制律，如PID控制、最优控制等在实际生产***中有着大量的应用。闭环控制律的引入，使得输入输出关系、***变量之间关系发生改变，并且闭环控制规律使得***对于外部干扰更具鲁棒性，使得故障处于早期阶段或幅值较小时，所带来的影响在***中被淹没，无法通过残差信号及时诊断，造成误报和漏报。这些限制了之前研究的诊断方法在实际中得应用，因此有必要在闭环前提下研究***状态监测和故障检测方法，以解决实际的生产安全性问题。目前基于数据的闭环***监测与故障检测研究成果十分有限，还处于初步的探索阶段，有待进一步深入研究。
本发明针对闭环***的过程监测与故障检测，并从工业***产品质量的角度出发，提出了一种结合慢特征分析和典型相关性分析的闭环***质量相关精细运行状态监测的算法。本方法通过提取过程数据与产品质量变量的相关性，实现对过程数据空间进行充分分解，分成与质量密切相关和与质量无关的两个空间，并建立模型对两个空间的动态和静态信息进行监测，该方法通过对动态和静态信息的分别监测可以实现对闭环***工况变化和真正故障发生的区分，并且通过质量相关的空间分解可以识别出异常具体发生在哪一个空间是否影响产品质量，提高了对***状态监测的精准，到目前为止，尚未见到与本发明相关的研究报道。
发明内容
本发明的目的在于针对工业闭环***过程故障检测，并且关注过程产品质量，提供了一种基于质量相关慢特征充分分解的闭环***精细运行状态监测方法。
本发明的目的是通过以下技术方案实现的：一种基于质量相关慢特征充分分解的闭环***精细运行状态监测方法，包含以下步骤：
(1)选择***正常过程测量数据与所关注的产品质量数据：对于一个工业闭环***，设其运行过程包含m个可测过程变量，在t时刻采样可以得到一个1×m 的向量x(t)＝(x1(t),x2(t),…,xm(t))，经过N次采样后得到正常过程下过程测量变量的数据矩阵XM＝(x(t),x(t+1),…,x(t+N))T，根据工业过程选择k个所关注的产品质量变量，经N次采样后，得到质量变量的数据矩阵Yk(∈N*k)。
(2)数据标准化：分别对过程数据矩阵XM和质量变量数据矩阵Yk按列减去该列的均值，并除以该列标准差进行标准化。
(3)对(2)中标准化后的过程数据矩阵X进行慢特征分析建模，获得初始慢特征。
这里考虑线性慢特征分析，每一个慢特征都可以看作是过程变量的线性组合，所以从过程数据矩阵X到慢特征S(t)＝[s1(t),s2(t),…,sm(t)]的映射表示成：
S(t)＝W×X (1)
其中W＝[w1,w2,…,wm]T是SFA需要优化的系数矩阵，m为慢特征的个数，与可测过程变量相同。求解慢特征的问题是使得所获得慢特征变化速率Δ(si)最小，求解慢特征问题转化成求解下列广义特征值问题：
AW＝BWΩ (2)
其中A表示输入的过程数据矩阵X的一阶差分矩阵的协方差矩阵，表示成B表示X矩阵的协方差矩阵<XXT>t，对应的，W＝[w1,w2,…,wm]T为特征向量构成的特征矩阵，m过程输入数据的变量个数；Ω为对应广义特征值构成的对角矩阵。利用两步奇异值分解(SVD)来求解慢特征，具体过程如下：
(3.1)首先对矩阵X进行白化，即对B矩阵的奇异值分解。矩阵X的协方差矩阵可表示成<XXT>t，对其进行奇异值分解(SVD)得：
<XXT>t＝UΛUT (3)
其中U是协方差矩阵<XXT>t的特征向量组成的矩阵，Λ是一个对角阵，其每一
个对角线上的元素就是一个特征值，则白化后得数据矩阵可表示成：
Z＝Λ-1/2UTX＝QX (4)
白化后数据矩阵Z得协方差满足<ZZT>t＝Q<XXT>tQT。
(3.2)求解线性SFA问题的目标等同于寻找一个矩阵P＝WQ-1使得S＝P*Z，并且使得S满足<SST>t＝P<ZZT>tPT＝I。
矩阵P可以通过对Z的差分矩阵的协方差矩阵进行奇异值分解求得：
由此可以计算得到初始慢特征系数矩阵为：
W＝PΛ-1/2UT (6)
初始慢特征为：
S＝WX＝PΛ-1/2UTX (7)
(4)利用(7)中获得的初始慢特征S与质量变量Yk进行典型相关性分析(CCA) 获取典型变量，具体过程如下：
对于初始慢特征S与质量变量Yk，定义u，v两个变量，其中u，v是始慢特征S与质量变量Yk的线性组合如下：
其中Ψ1＝[a1,a2,…,am]T，Ψ2＝[b1,b2,…,bk]T。求解典型变量u，v就是求解系数矩阵Ψ1，Ψ2使得u，v间的Pearson系数最大。转化成如下优化问题求解：
其中corr(u,v)表示u，v之间的Pearson系数，cov(u,v)为u，v的协方差，Var(u) 和Var(v)分别为u，v的方差。此问题通过构造Lagrangian等式来求解，求解对应最大Pearson系数的两个典型变量u，v。
同理，在此基础上继续求解寻找对应第二大Pearson系数的典型变量u2和v2：
一共可以求取n对典型变量对，n＝min(m,k)，对应的Pearson系数依次减小，即表示典型变量间的相关性越来越小。根据所求得的Pearson系数选取前p个对应 Pearson系数较大的典型相关变量对，舍去剩下对应Pearson系数较小的(n-p) 个典型相关变量对。利用前p个典型相关变量对，利用所求得关于慢特征S的线性组合的典型变量ui(i＝1,2,…,p)，组合成新的矩阵Uy＝[u1,u2,…,up]，其中 ui(i＝1,2,…,p)对应的Pearson系数依次减小，并且ui之间相互正交。
(5)将(4)中所获得的矩阵Uy利用偏最小二乘公式重构过程数据X空间，重构公式如下：
Eo为重构残差，即与Uy无关的部分。由此可将X空间分成与Yk相关的部分Uy和与 Yk无关的部分Xo两个部分，实现对X空间的质量充分分解。
X＝Uy+Xo (12)
(6)重新对Xo进行慢特征分析，得到慢特征，并按变化速度从小到大的顺序排列：
So＝WoXo＝[So1,So2,…,Som] (13)
Xo是原始过程数据减去质量变量Yk相关部分数据所获得的，所以获得的慢特征不是全部都有意义，选取满足的g个慢特征，E{}为求期望运算符，I是单位矩阵；
根据所获得慢特征的变化快慢程度，q个变化速度较慢特征组成矩阵 Sod＝[so1,so2,…,soq]，剩下(g-q)个变化速度较快的特征组成另一矩阵： Soe＝[so(q+1),so(q+2),…,sog]。
(7)建立静态和动态指标体系。
考虑到对于实际工业***,所关注的产品质量变量可能为一个或者两个以上,对于两种情况,建立不同的指标体系对***进行监测.
(7.1)对于单个质量变量的情况
由于质量变量只有一个,所以所提取到的典型变量也只有一个,所以对于与质量相关的空间直接利用典型变量Uy进行监测。其中静态指标为：
动态指标为：
对于动静态指标均采用shewhart控制图的方法进行监测，对于静态指标：
其中UTH1为的控制上限，LTH1为的控制下限，μ1为的均值，为其标准差，b1为控制图门限宽度。对于动态指标：
其中UTH2为的控制上限，LTH2为的控制下限，μ2为的均值，为其标准差，b2为控制图门限宽度。
对于与质量变量无关的空间，静态指标为：
动态指标为：
表示Sod的协方差矩阵特征值所构成的对角阵的逆，表示Soe的协方差矩阵特征值所构成的对角阵的逆，指标的控制限均通过核密度函数在0.95的置信水平下获得，分别为 分别为协方差矩阵的特征值构成的对角阵的逆，两个动态指标所对应的控制限均通过核密度函数在0.95的置信水平下计算得到，分别为
(7.2)对于有两个以上质量变量的情况
所获得的典型变量并不唯一，所以需要计算统计量来对分Uy进行监测，静态指标计算如下:
其中表示Uy的协方差矩阵对应的特征值构成的对角阵的逆，指标的控制限通过核密度函数在0.95的置信水平下获得分别为对Uy进行一次差分得到动态指标为：
为协方差矩阵的特征值构成的对角阵的逆，为其三个动态指标所对应的控制限，通过核密度函数在0.95的置信水平下计算得到。与质量变量无关的空间的指标计算方式与步骤7.1中单个质量变量的情况相同。
(8)利用动静态指标可以分别对闭环***的过程动态性和稳态操作情况进行监测。
对于与质量变量相关的空间:
1)静态指标超限,动态指标超限后又恢复到控制限以下,表明与质量变量相关的空间经过控制作用调节到达一个新的工况,质量变量出现短暂波动后又恢复到设定水平；
2)静态指标超限,动态指标也一直超限,表明与质量变量相关的空间出现故障,控制器调节失败,产品质量遭到破坏；
3)静态指标动态指标均未超限,则表明与质量变量相关的空间处于良好的运行状态,产品质量维持在设定值附近。
对于与质量变量无关的空间:
1)质量无关的静态指标与动态指标均短暂超限后回到控制限以下,表明在与质量变量无关的空间出现短暂异常,后通过闭环控制作用的调节,该空间又恢复到与异常发生前相同的稳态工况下运行；
2)与质量无关空间的静态指标一直超限,动态指标短暂超限后回到控制限以下,表明与质量变量无关的空间出现短暂异常,但是通过闭环控制作用的调节该空间恢复到一个异常发生前并不一样的稳态工况下运行,静态性能已经发生改变,但是动态性能仍维持正常；
3)质量无关的静态指标与动态指标均超限,表明该空间出现故障,且该空间相关控制器调节失败,静动态性能都遭到破坏；
4)质量无关的静态指标与动态指标均维持正常,表明该空间保持一个良好的运行状态。
本发明的有益效果为：本发明可以将工业闭环***划分成与产品质量相关和与产品质量无关的两个空间，不仅可以有效区分故障具体发生在哪一个空间，进而判断故障是否影响产品质量，而且通过对***动静态信息的监测，可以有效地识别工况的变化与真正故障的区分，实现闭环***的精细运行状态监测。
附图说明
图1是本发明的流程图；
图2是本方法具体应用的闭环TE过程工艺流程图；
图3为TE过程提供的故障类展示图；
图4是本方法应用在TE过程故障1情况下的静态指标监测结果；
图5是本方法应用在TE过程故障1情况下的动态指标监测结果；
图6是TE过程故障1情况下对应的产品质量A与正常情况下产品质量A对比图；
图7是TE过程故障1情况下对应的产品质量B与正常情况下产品质量B对比图；
图8是TE过程故障1情况下对应的产品质量C与正常情况下产品质量C对比图；
图9是本方法应用在TE过程故障4情况下的静态指标监测结果；
图10是本方法应用在TE过程故障4情况下的动态指标监测结果；
图11是TE过程故障4情况下对应的产品质量A与正常情况下产品质量A 对比图；
图12是TE过程故障4情况下对应的产品质量A与正常情况下产品质量A 对比图；
图13是TE过程故障4情况下对应的产品质量A与正常情况下产品质量A 对比图。
具体实施方式
下面结合附图及具体实例，对本发明作进一步详细说明。
Tenessee Eastman(TE)过程是由Downs等人基于Tenessee Eastman化学公司某实际化工生产过程提出的一个仿真***,在过程***工程领域的研究中,TE 过程是一个常用的标准问题(Benchmark problem),其较好地模拟了实际复杂工业过程***的许多典型特征,因此被作为仿真例子广泛应用于控制、优化、过程监控与故障诊断的研究中。本研究中，采用的是Lyman和Georgakis提出的全厂控制策略的闭环TE过程。TE过程包含四种气体原料A、C、D和E,两种液态产物G和H,还包含副产物F和惰性气体B。
TE过程包括五个主要单元：反应器、冷凝器、压缩机、分离器和汽提塔，如图2所示，共包括41个测量变量和12个控制变量。如图3，TE过程提供了 21个故障类,可供使用。
如图1所示，本发明包括以下步骤：
(1)数据选择:选择TE过程正常运行时的过程测量数据XMEAS(1-22)和 XMV(1-11)共33个变量作为过程输入变量X，XMEAS(29-31)三个变量所关注的产品量变量Y，每次采样可以得到一个向量x(t)＝(x1(t),x2(t),…,xm(t))维度为1×m 的，经过N次采样后得到过程测量变量X＝(x(t),x(t+1),…,x(t+N))T，产品质量变量，经N次采样后，得到质量变量的数据矩阵Yk(∈N*k)。
(2)数据标准化：分别对过程数据矩阵X和质量变量数据矩阵Yk按列去均值并除以标准差进行标准化。
(3)对(2)中标准化后的过程数据矩阵X进行慢特征分析建模，获得初始慢特征。
这里考虑线性慢特征分析，每一个慢特征都可以看作是过程变量的线性组合，所以从过程数据矩阵到慢特征S(t)＝[s1(t),s2(t),…,sm(t)]的映射可以表示成：
S(t)＝W×X
其中W＝[w1,w2,…,wm]T是SFA需要优化的系数矩阵，m为慢特征的个数，与过程变量数相同。求解慢特征的问题是使得所获得慢特征变化速率Δ(si)最小，求解慢特征问题可以转化成求解下列广义特征值问题：
AW＝BWΩ
其中A表示输入数据X的一阶差分矩阵的协方差矩阵，可表示成B 表示X矩阵的协方差矩阵<XXT>t，W＝[w1,w2,…,wm]T为特征向量构成的特征矩阵，m为过程输入数据的变量个数，Ω为对应广义特征值构成的对角矩阵，同样也是。同样可以利用两步奇异值分解(SVD)来求解慢特征，具体过程如下：
(3.1)首先对矩阵X进行白化，即对B矩阵的奇异值分解。矩阵X的协方差矩阵可表示成<XXT>t，对其进行奇异值分解(SVD)得：
<XXT>t＝UΛUT
其中U是协方差矩阵<XXT>t的特征向量组成的矩阵，Λ是一个对角阵，其每一
个对角线上的元素就是一个特征值，则白化后得数据矩阵可表示成：
Z＝Λ-1/2UTX＝QX
白化后数据矩阵Z得协方差满足<ZZT>t＝Q<XXT>tQT。
(4.2)求解线性SFA问题的目标等同于寻找一个矩阵P＝WQ-1使得S＝P*Z，并且使得S满足<SST>t＝P<ZZT>tPT＝I。
矩阵P可以通过对Z的差分矩阵的协方差矩阵进行奇异值分解求得：
由此可以计算得到初始慢特征系数矩阵为：
W＝PΛ-1/2UT
初始慢特征为：
S＝WX＝PΛ-1/2UTX
(5)利用步骤(4)中获得的初始慢特征S与质量变量Yk进行典型相关性分析 (CCA)获取典型变量，具体过程如下：
(5.1)对于初始慢特征S与质量变量Yk，定义u，v两个变量，其中u，v 是始慢特征S与质量变量Yk的线性组合如下：
求解典型变量u，v就是求解系数矩阵A，B使得u，v间的Pearson系数最大。可转化成如下优化问题求解：
Maximizeu2
其中corr(u,v)表示Pearson系数，cov(u,v)为u，v的协方差，Var(u)和Var(v)分别u，v的方差。此问题可以通过构造Lagrangian等式求解，求解对应最大Pearson 系数的两个典型变量u，v。
(5.1)同理，可以在此基础上继续求解寻找对应第二大Pearson系数的典型变量u2和v2：
一共可以求取n对典型变量对，n＝min(m,k)，对应的Pearson系数依次减小，即典型变量间的相关性越来越小。根据需要一般根据所求得的Pearson系数选取前 p个对应最大相关系数的典型相关变量对，舍去剩下对应Pearson系数很小的(n-p) 个典型相关变量对。利用前p个典型相关变量对，利用对应慢特征S的线性组合的典型变量组合成新的矩阵Uy＝[u1,u2,…,up]，ui(i＝1,2,…,p)为求得的典型变量，对应的Pearson系数依次减小，并且ui之间互不相关。
(6)将(5)中所获得的典型变量Uy利用偏最小二乘公式重构回Xd空间，重构公式如下：
X＝UyPT+E
PT＝(Uy TUy)-1Uy TXd
并将X空间分成与Yk相关的部分Uy和与Yk无关的部分Xo两个部分，实现对X空间的质量充分分解。
X＝Uy+Xo
(7)重新对Xo进行慢特征分析，得到慢特征，并按变化速度从小到大的顺序排列：
So＝WoXo＝[So1,So2,…,Som]
Xo是原始过程数据减去质量变量Yk相关部分数据所获得的，所以获得的慢特征不是全部都有意义，选取满足的g个慢特征。并根据所获得慢特征的变化快慢程度，选取q个变化速度较慢特征组成Sod＝[so1,so2,…,soq]，剩下(g-q) 个变化速度较快的特征Soe＝[so(q+1),so(q+2),…,sog]。
(8)建立静态和动态指标体系。
考虑到对于实际工业***,所关注的产品质量变量可能为一个或者多个,对于两种情况,建立不同的指标体系对***进行监测.
(8.1)对于单个质量变量的情况
由于质量变量只有一个,所以所提取到的典型变量也只有一个,所以对于与质量相关的空间直接利用典型变量Uy进行监测。其中静态指标为：
动态指标为：
对于动静态指标均采用shewhart控制图的方法进行监测，以静态指标为例：
其中UTH为的控制上限，LTH为控制下限，μ为的均值，为其标准差，b为控制图门限宽度。
对于与质量变量无关的空间，静态指标为：
动态指标为：
表示Sod的协方差矩阵特征值所构成的对角阵的逆，表示Soe的协方差矩阵特征值所构成的对角阵的逆，指标的控制限均通过核密度函数在0.95的置信水平下获得分别为 分别为协方差矩阵的特征值构成的对角阵的逆，分别为两个动态指标所对应的控制限，通过核密度函数在0.95的置信水平下计算得到。
(8.2)对于有多个质量变量的情况
所获得的典型变量并不唯一，所以需要计算统计量来对分Uy进行监测，静态指标计算如下:
其中表示Uy的协方差矩阵对应的特征值构成的对角阵的逆，指标的控制限均通过核密度函数在0.95的置信水平下获得分别为对Uy进行一次差分得到动态指标为：
为协方差矩阵的特征值构成的对角阵的逆，为其三个动态指标所对应的控制限，通过核密度函数在0.95的置信水平下计算得到,
与质量变量无关的空间的指标计算方式同上。
利用动静态指标可以分别对闭环***的过程动态性和稳态操作情况进行监测。对于与质量变量相关的空间:
1)静态指标超限,动态指标超限后又恢复到控制限以下,表明与质量变量相关的空间经过控制作用调节到达一个新的工况,质量变量出现短暂波动后又恢复到设定水平；
2)静态指标超限,动态指标也一直超限,表明与质量变量相关的空间出现故障,控制器调节失败,产品质量遭到破坏；
3)静态指标动态指标均未超限,则表明与质量变量相关的空间处于良好的运行状态,产品质量维持再设定值附近。
对于与质量变量无关的空间:
1)质量无关的静态指标与动态指标均短暂超限后回到控制限以下,表明在与质量变量无关的空间出现短暂异常,后通过闭环控制作用的调节,该空间又恢复到与异常发生前相同的稳态工况下运行；
2)与质量无关空间的静态指标一直超限,动态指标短暂超限后回到控制限以下,表明与质量变量无关的空间出现短暂异常,但是通过闭环控制作用的调节该空间恢复到一个异常发生前并不一样的稳态工况下运行,静态性能已经发生改变,但是动态性能仍维持正常；
3)质量无关的静态指标与动态指标均超限,表明该空间出现故障,且该空间相关控制器调节失败,静动态性能都遭到破坏；
4)质量无关的静态指标与动态指标均维持正常,表明该空间保持一个良好的运行状态。
两个空间相互独立,相应的动静态指标分别指示各自空间的动静态信息的变化,可以帮助我们判断故障及异常具体发生在哪个空间,是否会影响产品质量。
图4、5本方法用于TE过程故障1数据的监测结果，图6、7、8是所选择的三个产品质量变量正常情况下变化和在故障1情况下的变化对比。根据图可以看出，与产品质量y相关的空间，静态指标和动态指标超限一段时间后均恢复到控制限以下，说明该空间受到故障的短暂干扰，通过控制器调节一段时间后静动态特性均恢复到初始水平，产品质量也恢复到初始设定值；与产品质量无关的空间，静态指标一直超限，但趋于平稳，动态指标超限一段时间后恢复到控制限以下，说明该空间受到干扰后通过调节恢复到一个新的稳态工况下。图9、10是本方法用于TE过程故障4数据的监测结果，从图中可以看出与产品质量y相关的空间静动态指标均未超限，说明该空间运行正常，产品质量未产生变化，从图11、12、13中各个产品质量正常情况下和故障四情况下的变化对比也可以看出监测结果的判断是正确的；在与产品质量无关的空间，静态指标一直超限，动态指标短暂超限后恢复到控制限以下，说明故障发生在与产品质量无关的闭环中，通过控制器的调节又达到一个新的稳态工况。传统的慢特征分析只能区分工况异常和故障发生，而通过本方法可以进一步区分故障所发生的空间以及是否影响产品质量。
Claims (1)
1.一种基于质量相关慢特征充分分解的闭环***精细运行状态监测方法，其特征在于，包含以下步骤：
(1)选择***正常过程测量数据与所关注的产品质量数据：对于一个工业闭环***，设其运行过程包含m个可测过程变量，在t时刻采样可以得到一个1×m的向量x(t)＝(x1(t),x2(t),…,xm(t))，经过N次采样后得到正常过程下过程测量变量的数据矩阵XM＝(x(t),x(t+1),…,x(t+N))T，根据工业过程选择k个所关注的产品质量变量，经N次采样后，得到质量变量的数据矩阵Yk(∈N*k)。
(2)数据标准化：分别对过程数据矩阵XM和质量变量数据矩阵Yk按列减去该列的均值，并除以该列标准差进行标准化。
(3)对(2)中标准化后的过程数据矩阵X进行慢特征分析建模，获得初始慢特征。
这里考虑线性慢特征分析，每一个慢特征都可以看作是过程变量的线性组合，所以从过程数据矩阵X到慢特征S(t)＝[s1(t),s2(t),…,sm(t)]的映射表示成：
S(t)＝W×X (1)
其中W＝[w1,w2,…,wm]T是SFA需要优化的系数矩阵，m为慢特征的个数，与可测过程变量相同。求解慢特征的问题是使得所获得慢特征变化速率Δ(si)最小，求解慢特征问题转化成求解下列广义特征值问题：
AW＝BWΩ (2)
其中A表示输入的过程数据矩阵X的一阶差分矩阵的协方差矩阵，表示成B表示X矩阵的协方差矩阵<XXT>t，对应的，W＝[w1,w2,…,wm]T为特征向量构成的特征矩阵，m过程输入数据的变量个数；Ω为对应广义特征值构成的对角矩阵。利用两步奇异值分解(SVD)来求解慢特征，具体过程如下：
(3.1)首先对矩阵X进行白化，即对B矩阵的奇异值分解。矩阵X的协方差矩阵可表示成<XXT>t，对其进行奇异值分解(SVD)得：
<XXT>t＝UΛUT (3)
其中U是协方差矩阵<XXT>t的特征向量组成的矩阵，Λ是一个对角阵，其每一个对角线上的元素就是一个特征值，则白化后得数据矩阵可表示成：
Z＝Λ-1/2UTX＝QX (4)
白化后数据矩阵Z得协方差满足<ZZT>t＝Q<XXT>tQT。
(3.2)求解线性SFA问题的目标等同于寻找一个矩阵P＝WQ-1使得S＝P*Z，并且使得S满足<SST>t＝P<ZZT>tPT＝I。
矩阵P可以通过对Z的差分矩阵的协方差矩阵进行奇异值分解求得：
由此可以计算得到初始慢特征系数矩阵为：
W＝PΛ-1/2UT (6)
初始慢特征为：
S＝WX＝PΛ-1/2UTX (7)
(4)利用(7)中获得的初始慢特征S与质量变量Yk进行典型相关性分析(CCA)获取典型变量，具体过程如下：
对于初始慢特征S与质量变量Yk，定义u，v两个变量，其中u，v是始慢特征S与质量变量Yk的线性组合如下：
其中Ψ1＝[a1,a2,…,am]T，Ψ2＝[b1,b2,…,bk]T。求解典型变量u，v就是求解系数矩阵Ψ1，Ψ2使得u，v间的Pearson系数最大。转化成如下优化问题求解：
Subjectto：
其中corr(u,v)表示u，v之间的Pearson系数，cov(u,v)为u，v的协方差，Var(u)和Var(v)分别为u，v的方差。此问题通过构造Lagrangian等式来求解，求解对应最大Pearson系数的两个典型变量u，v。
同理，在此基础上继续求解寻找对应第二大Pearson系数的典型变量u2和v2：
一共可以求取n对典型变量对，n＝min(m,k)，对应的Pearson系数依次减小，即表示典型变量间的相关性越来越小。根据所求得的Pearson系数选取前p个对应Pearson系数较大的典型相关变量对，舍去剩下对应Pearson系数较小的(n-p)个典型相关变量对。利用前p个典型相关变量对，利用所求得关于慢特征S的线性组合的典型变量ui(i＝1,2,…,p)，组合成新的矩阵Uy＝[u1,u2,…,up]，其中ui(i＝1,2,…,p)对应的Pearson系数依次减小，并且ui之间相互正交。
(5)将(4)中所获得的矩阵Uy利用偏最小二乘公式重构过程数据X空间，重构公式如下：
Eo为重构残差，即与Uy无关的部分。由此可将X空间分成与Yk相关的部分Uy和与Yk无关的部分Xo两个部分，实现对X空间的质量充分分解。
X＝Uy+Xo (12)
(6)重新对Xo进行慢特征分析，得到慢特征，并按变化速度从小到大的顺序排列：
So＝WoXo＝[So1,So2,…,Som] (13)
Xo是原始过程数据减去质量变量Yk相关部分数据所获得的，所以获得的慢特征不是全部都有意义，选取满足的g个慢特征，E{}为求期望运算符，I是单位矩阵；
根据所获得慢特征的变化快慢程度，q个变化速度较慢特征组成矩阵Sod＝[so1,so2,…,soq]，剩下(g-q)个变化速度较快的特征组成另一矩阵：Soe＝[so(q+1),so(q+2),…,sog]。
(7)建立静态和动态指标体系。
考虑到对于实际工业***,所关注的产品质量变量可能为一个或者两个以上,对于两种情况,建立不同的指标体系对***进行监测.
(7.1)对于单个质量变量的情况
由于质量变量只有一个,所以所提取到的典型变量也只有一个,所以对于与质量相关的空间直接利用典型变量Uy进行监测。其中静态指标为：
动态指标为：
对于动静态指标均采用shewhart控制图的方法进行监测，对于静态指标：
其中UTH1为的控制上限，LTH1为的控制下限，μ1为的均值，为其标准差，b1为控制图门限宽度。对于动态指标：
其中UTH2为的控制上限，LTH2为的控制下限，μ2为的均值，为其标准差，b2为控制图门限宽度。
对于与质量变量无关的空间，静态指标为：
动态指标为：
表示Sod的协方差矩阵特征值所构成的对角阵的逆，表示Soe的协方差矩阵特征值所构成的对角阵的逆，指标的控制限均通过核密度函数在0.95的置信水平下获得，分别为 分别为协方差矩阵的特征值构成的对角阵的逆，两个动态指标所对应的控制限均通过核密度函数在0.95的置信水平下计算得到，分别为
(7.2)对于有两个以上质量变量的情况
所获得的典型变量并不唯一，所以需要计算统计量来对分Uy进行监测，静态指标计算如下:
其中表示Uy的协方差矩阵对应的特征值构成的对角阵的逆，指标的控制限通过核密度函数在0.95的置信水平下获得分别为对Uy进行一次差分得到动态指标为：
为协方差矩阵的特征值构成的对角阵的逆，为其三个动态指标所对应的控制限，通过核密度函数在0.95的置信水平下计算得到。与质量变量无关的空间的指标计算方式与步骤7.1中单个质量变量的情况相同。
(8)利用动静态指标可以分别对闭环***的过程动态性和稳态操作情况进行监测。
对于与质量变量相关的空间:
1)静态指标超限,动态指标超限后又恢复到控制限以下,表明与质量变量相关的空间经过控制作用调节到达一个新的工况,质量变量出现短暂波动后又恢复到设定水平；
2)静态指标超限,动态指标也一直超限,表明与质量变量相关的空间出现故障,控制器调节失败,产品质量遭到破坏；
3)静态指标动态指标均未超限,则表明与质量变量相关的空间处于良好的运行状态,产品质量维持在设定值附近。
对于与质量变量无关的空间:
1)质量无关的静态指标与动态指标均短暂超限后回到控制限以下,表明在与质量变量无关的空间出现短暂异常,后通过闭环控制作用的调节,该空间又恢复到与异常发生前相同的稳态工况下运行；
2)与质量无关空间的静态指标一直超限,动态指标短暂超限后回到控制限以下,表明与质量变量无关的空间出现短暂异常,但是通过闭环控制作用的调节该空间恢复到一个异常发生前并不一样的稳态工况下运行,静态性能已经发生改变,但是动态性能仍维持正常；
3)质量无关的静态指标与动态指标均超限,表明该空间出现故障,且该空间相关控制器调节失败,静动态性能都遭到破坏；
4)质量无关的静态指标与动态指标均维持正常,表明该空间保持一个良好的运行状态。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810771664.6A CN109143995B (zh) | 2018-07-13 | 2018-07-13 | 质量相关慢特征分解的闭环***精细运行状态监测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810771664.6A CN109143995B (zh) | 2018-07-13 | 2018-07-13 | 质量相关慢特征分解的闭环***精细运行状态监测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109143995A true CN109143995A (zh) | 2019-01-04 |
| CN109143995B CN109143995B (zh) | 2020-09-01 |
Family
ID=64800694
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810771664.6A Active CN109143995B (zh) | 2018-07-13 | 2018-07-13 | 质量相关慢特征分解的闭环***精细运行状态监测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109143995B (zh) |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110222393A (zh) * | 2019-05-28 | 2019-09-10 | 浙江大学 | 基于细粒度风电发电状态划分的风机叶片结冰异常监测方法 |
| CN111401940A (zh) * | 2020-03-05 | 2020-07-10 | 杭州网易再顾科技有限公司 | 特征预测方法、装置、电子设备及存储介质 |
| CN111413949A (zh) * | 2020-03-30 | 2020-07-14 | 南京富岛信息工程有限公司 | 一种降低工业过程故障预警误报率的方法 |
| CN111624979A (zh) * | 2020-05-18 | 2020-09-04 | 浙江大学 | 一种基于慢特征分析的工业闭环控制回路多振荡检测与溯源方法 |
| CN111949003A (zh) * | 2020-07-17 | 2020-11-17 | 浙江浙能技术研究院有限公司 | 一种基于SFA与Hellinger距离的闭环控制回路性能评价方法 |
| CN112199409A (zh) * | 2020-08-17 | 2021-01-08 | 浙江中控技术股份有限公司 | 一种用于催化重整装置实时工况的监测方法及装置 |
| CN112541017A (zh) * | 2020-12-02 | 2021-03-23 | 中海油信息科技有限公司 | 一种工业生产过程状态监测方法、装置、设备及存储介质 |
| CN117048113A (zh) * | 2023-10-13 | 2023-11-14 | 江苏爱箔乐铝箔制品有限公司 | 铝箔餐盒冲压的动态精度检测方法及*** |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104077608A (zh) * | 2014-06-11 | 2014-10-01 | 华南理工大学 | 基于稀疏编码慢特征函数的行为识别方法 |
| CN104537260A (zh) * | 2015-01-14 | 2015-04-22 | 清华大学 | 基于缓慢特征回归的动态软测量方法和*** |
| CN105678723A (zh) * | 2015-12-29 | 2016-06-15 | 内蒙古科技大学 | 基于稀疏分解和差分图像的多聚焦图像融合方法 |
-
2018
- 2018-07-13 CN CN201810771664.6A patent/CN109143995B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104077608A (zh) * | 2014-06-11 | 2014-10-01 | 华南理工大学 | 基于稀疏编码慢特征函数的行为识别方法 |
| CN104537260A (zh) * | 2015-01-14 | 2015-04-22 | 清华大学 | 基于缓慢特征回归的动态软测量方法和*** |
| CN105678723A (zh) * | 2015-12-29 | 2016-06-15 | 内蒙古科技大学 | 基于稀疏分解和差分图像的多聚焦图像融合方法 |
Non-Patent Citations (2)
| Title |
|---|
| SHANG, C.等: "Concurrent monitoring of operating condition deviations and process dynamics anomalies with slow feature analysis", 《AICHE JOURNAL》 * |
| 蒋昕祎 等: "基于核慢特征回归与互信息的常压塔软测量建模", 《化工学报》 * |
Cited By (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110222393A (zh) * | 2019-05-28 | 2019-09-10 | 浙江大学 | 基于细粒度风电发电状态划分的风机叶片结冰异常监测方法 |
| CN111401940A (zh) * | 2020-03-05 | 2020-07-10 | 杭州网易再顾科技有限公司 | 特征预测方法、装置、电子设备及存储介质 |
| CN111413949B (zh) * | 2020-03-30 | 2023-12-01 | 南京富岛信息工程有限公司 | 一种降低工业过程故障预警误报率的方法 |
| CN111413949A (zh) * | 2020-03-30 | 2020-07-14 | 南京富岛信息工程有限公司 | 一种降低工业过程故障预警误报率的方法 |
| CN111624979A (zh) * | 2020-05-18 | 2020-09-04 | 浙江大学 | 一种基于慢特征分析的工业闭环控制回路多振荡检测与溯源方法 |
| CN111624979B (zh) * | 2020-05-18 | 2021-07-06 | 浙江大学 | 一种基于慢特征分析的工业闭环控制回路多振荡检测与溯源方法 |
| CN111949003A (zh) * | 2020-07-17 | 2020-11-17 | 浙江浙能技术研究院有限公司 | 一种基于SFA与Hellinger距离的闭环控制回路性能评价方法 |
| CN111949003B (zh) * | 2020-07-17 | 2021-09-03 | 浙江浙能技术研究院有限公司 | 一种基于SFA与Hellinger距离的闭环控制回路性能评价方法 |
| CN112199409A (zh) * | 2020-08-17 | 2021-01-08 | 浙江中控技术股份有限公司 | 一种用于催化重整装置实时工况的监测方法及装置 |
| CN112199409B (zh) * | 2020-08-17 | 2024-06-11 | 浙江中控技术股份有限公司 | 一种用于催化重整装置实时工况的监测方法及装置 |
| CN112541017A (zh) * | 2020-12-02 | 2021-03-23 | 中海油信息科技有限公司 | 一种工业生产过程状态监测方法、装置、设备及存储介质 |
| CN117048113B (zh) * | 2023-10-13 | 2023-12-22 | 江苏爱箔乐铝箔制品有限公司 | 铝箔餐盒冲压的动态精度检测方法及*** |
| CN117048113A (zh) * | 2023-10-13 | 2023-11-14 | 江苏爱箔乐铝箔制品有限公司 | 铝箔餐盒冲压的动态精度检测方法及*** |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109143995B (zh) | 2020-09-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109143995A (zh) | 一种基于质量相关慢特征充分分解的闭环***精细运行状态监测方法 | |
| Shang et al. | Fault detection based on augmented kernel Mahalanobis distance for nonlinear dynamic processes | |
| Joe Qin | Statistical process monitoring: basics and beyond | |
| CA2679128C (en) | Method and system of using inferential measurements for abnormal event detection in continuous industrial processes | |
| Peng et al. | Quality‐related process monitoring based on total kernel PLS model and its industrial application | |
| Chiang et al. | Fault detection and diagnosis in industrial systems | |
| Ait-Izem et al. | On the application of interval PCA to process monitoring: A robust strategy for sensor FDI with new efficient control statistics | |
| Van den Kerkhof et al. | Dynamic model-based fault diagnosis for (bio) chemical batch processes | |
| He et al. | A novel process monitoring and fault detection approach based on statistics locality preserving projections | |
| CN108664002B (zh) | 一种面向质量的非线性动态过程监控方法 | |
| JP4046309B2 (ja) | プラント監視装置 | |
| Jiang et al. | Weighted kernel principal component analysis based on probability density estimation and moving window and its application in nonlinear chemical process monitoring | |
| Rotem et al. | Ethylene compressor monitoring using model‐based PCA | |
| Wang et al. | Fault isolation based on residual evaluation and contribution analysis | |
| Wang et al. | Data-driven sensor fault diagnosis systems for linear feedback control loops | |
| Angeli | Online expert systems for fault diagnosis in technical processes | |
| Ralston et al. | Computer-based monitoring and fault diagnosis: a chemical process case study | |
| CN109145256B (zh) | 一种基于规范变量非线性主成分分析的过程监测方法 | |
| Slišković et al. | Multivariate statistical process monitoring | |
| Cheng et al. | Rebooting kernel CCA method for nonlinear quality-relevant fault detection in process industries | |
| CN111796576B (zh) | 一种基于双核t分布随机近邻嵌入的过程监测可视化方法 | |
| Xiong et al. | Multivariate statistical process monitoring of an industrial polypropylene catalyzer reactor with component analysis and kernel density estimation | |
| Hallgrímsson et al. | Improved process diagnosis using fault contribution plots from sparse autoencoders | |
| Linzhe et al. | A nonlinear quality-relevant process monitoring method with kernel input-output canonical variate analysis | |
| CN109283912B (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 |