基于GSFA-GNPE的动态-静态联合指标间歇过程监控

doi:10.16183/j.cnki.jsjtu.2020.290

基于GSFA-GNPE的动态-静态联合指标间歇过程监控

赵小强^,^a^,^b^,^c, 牟淼^a

a.电气工程与信息工程学院, 兰州理工大学, 兰州 730050

b.甘肃省工业过程先进控制重点实验室, 兰州理工大学, 兰州 730050

c.电气与控制工程国家级实验教学示范中心, 兰州理工大学, 兰州 730050

Batch Process Monitoring with Dynamic-Static Joint Indicator Based on GSFA-GNPE

ZHAO Xiaoqiang^,^a^,^b^,^c, MOU Miao^a

a. College of Electrical and Information Engineering, Lanzhou University of Technology, Lanzhou 730050, China

b. Key Laboratory of Advanced Control for Industrial Processes of Gansu Province, Lanzhou University of Technology, Lanzhou 730050, China

c. National Demonstration Center for Experimental Electrical and Control Engineering Education, Lanzhou University of Technology, Lanzhou 730050, China

责任编辑: 石易文

收稿日期: 2020-09-14

基金资助:

国家自然科学基金(61763029)
国防基础科研(JCKY2018427C002)
甘肃省高等学校产业支撑引导(2019C-05)
甘肃省工业过程先进控制重点实验室开放基金资助项目(2019KFJJ01)

Received: 2020-09-14

作者简介 About authors

赵小强(1969-),男,陕西省宝鸡市人,教授,博士生导师,研究方向为过程监控和故障诊断、生产调度及数据挖掘.E-mail:xqzhao@lut.edu.cn.

摘要

传统的过程监控方法忽略了变量间的时序相关性,且没有区分变量间的动态关系与静态关系,从而导致监控效果不佳.针对此问题,本文提出一种基于全局慢特征分析(GSFA)-全局邻域保持嵌入(GNPE)的动态-静态联合指标间歇过程监控方法,该方法可以有效提取动态全局特征和静态全局特征.首先,对过程变量的动态特性和静态特性进行评估,把自相关和互相关性较弱的变量视为静态变量,剩余变量视为动态变量;其次,分别对动态子空间和静态子空间构建GSFA和GNPE模型;然后,对来自每个子空间的统计信息使用贝叶斯推理进行组合,以得出混合模型的联合指标实现过程监控;最后,将所提算法应用于数值算例和青霉素发酵仿真过程进行仿真验证.结果表明,GSFA-GNPE算法相较于其他算法的故障检测效果更好.

关键词： 间歇过程; 过程监控; 慢特征分析; 邻域保持嵌入; 全局-局部; 贝叶斯推断

Abstract

Traditional process monitoring methods ignore the time-series correlation between variables, and do not distinguish the dynamic relationship and static relationship between variables, resulting in poor monitoring effect. To solve these problems, a dynamic-static joint indicator monitoring method of batch process based on global slow feature analysis(GSFA)-global neighborhood preserving embedding (GNPE) is proposed in this paper, which can effectively extract dynamic global features and static global features. First, the dynamic and static characteristics of the process variables are evaluated. Variables with weak autocorrelation and cross-correlation are regarded as static variables, and the remaining variables are regarded as dynamic ones. Next, the GSFA and GNPE models are constructed for dynamic and static subspaces, respectively. Finally, the statistical information from each subspace is combined by using Bayesian inference to obtain the joint indicator of the mixed model to realize process monitoring. Finally, the proposed algorithm is applied to a numerical example and the penicillin fermentation simulation process for simulation verification. The results show that the proposed GSFA-GNPE algorithm has better fault detection effects than other algorithms.

Keywords： batch process; process monitoring; slow feature analysis; neighborhood preserving embedding; globality-locality; Bayesian inference

PDF (2334KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

赵小强, 牟淼. 基于GSFA-GNPE的动态-静态联合指标间歇过程监控[J]. 上海交通大学学报（自然版）, 2021, 55(11): 1417-1428 doi:10.16183/j.cnki.jsjtu.2020.290

ZHAO Xiaoqiang, MOU Miao. Batch Process Monitoring with Dynamic-Static Joint Indicator Based on GSFA-GNPE[J]. Journal of shanghai Jiaotong University, 2021, 55(11): 1417-1428 doi:10.16183/j.cnki.jsjtu.2020.290

随着对生产过程安全性和产品质量要求的日益增长,过程监控受到越来越多的关注^[1,2].近年来,多元统计过程监控(MSPM)方法被广泛应用于过程监控,其主要思想是将高维数据投影到低维空间,并保留原始数据的主要信息^[3,4,5].

传统的MSPM方法, 如主成分分析(PCA)和邻域保持嵌入(NPE)等^[6,7,8,9], 假定在不同时间的样本是独立的,即过程当前时刻的状态不受之前时刻的影响,然而在现代工业中这种假设是不成立的,因此传统的监控模型无法准确表达过程数据的动态特性.为解决这一问题,文献[10]利用包含当前和过去时刻采样值的增广矩阵进行建模,提出了动态主元分析(DPCA)方法.然而,DPCA方法本质是一种全局结构保持方法,忽略了局部信息的提取,为此,NPE等局部流形算法被用于过程监控.NPE通过保持数据的局部邻域结构来实现数据降维,其假设数据样本满足独立分布,主要关注数据的空间结构特征,能有效提取数据的静态信息,但NPE算法忽视了数据随时间变化的特征信息,可能会造成检测效果不佳.为了有效提取过程数据的动态信息和局部信息,文献[11]在邻域保持嵌入算法保持数据局部结构的基础上引入时序扩展,提出一种基于时序扩展的邻域保持嵌入(TNPE)算法用于动态过程建模,一定程度上提高了动态过程的故障检测效果.但是,上述方法假设生产过程的所有变量均具有动态特性,并没有综合考虑生产过程中的静态变量.同时,NPE作为流形算法,其降维时主要考虑数据的局部结构特征会造成全局信息的丢失.因此,有必要为NPE算法加入全局目标函数加以约束,使其能充分提取过程信息.

近年来,慢特征分析(SFA)被用于过程监控^{[12,13,14,15]},SFA的主要思想是从时间序列中提取缓慢变化的特征^[16],能够有效提取表征过程数据动态变化的慢特征.文献[17]将SFA方法扩展到质量相关的故障检测领域,提出了一种基于SFA的工业过程质量相关故障检测算法.文献[18]将核方法应用于慢特征分析,提出了一种核慢特征判别分析(KSFDA)和支持向量数据描述(SVDD)的故障检测方法,对于非线性过程具有较好的检测效果.然而在实际工业过程中,一些过程变量可能是动态的,而另外一些则是静态的.尽管基于SFA的故障检测算法在动态工程监控中取得了长足的进步,但SFA未能同时考虑过程的动态特性和静态特性.对于动态变量,动态特征能够通过变量的动态变化准确地进行提取;而对于静态变量,动态信息则是由随机噪声引起的,提取的动态特征是与噪声相关的特征,与过程数据无关.此外,随机噪声可能会随着基准数据集的不同而变化,因此不需要在静态变量中进行动态特征提取.同时,将SFA应用于间歇过程监控时,由于间歇过程在进行数据展开时会破坏其原始时间结构信息,需要根据最近邻准则构造伪时间序列以提供时间结构信息,这就使得SFA等价于一种局部近邻数据关系保持算法,而忽略了全局结构信息.因而,需要将全局结构分析融入SFA,使其能充分提取过程信息.

将所有过程变量都视为静态变量用静态方法,或都视为动态变量用动态方法进行监控,都会导致故障监控效果不佳,且单一考虑全局信息或局部信息都会造成过程信息丢失.鉴于SFA和NPE在处理动态数据和静态数据方面各自的优势,提出一种基于全局慢特征分析(GSFA)-全局邻域保持嵌入(GNPE)的动态-静态联合指标间歇过程监控方法.首先,计算每个变量的自相关系数和互相关系数,分别求得当前时间和过去时间的变量之间的相关性以及当前时间的变量和过去的另一个变量之间的相关性.在此基础上,将过程变量分为动态子空间和静态子空间.在动态子空间中建立GSFA模型,提取动态全局信息;在静态子空间中建立GNPE模型,提取静态全局信息.最后,将两个子空间得出的监视统计信息进行贝叶斯推断,得出联合指标以进行故障检测.本文方法通过将过程变量区分为动态变量和静态变量,在不同的变量空间中采用不同的方法进行建模,可以充分提取过程信息,有效提高对于间歇过程的监控效果.

1 基础算法

1.1 慢特征分析算法

给定一个在时间t∈ $[t_{0}, t_{1}]$ 上连续变化的m维输入x(t)= ${[x_{1} (t) x_{2} (t) \dots x_{m} (t)]}^{T}$ ,SFA的目标函数是找到一个映射函数g(x),使得其输出信号s_i(t)=g_i(x(t))(i=1,2,…,m)具有最慢的时间变化,其中g(x)= ${[g_{1} (x) g_{2} (x) \dots g_{m} (x)]}^{T}$ 和s(t)= ${[s_{1} (t) s_{2} (t) \dots s_{m} (t)]}^{T}$ 分别为映射函数和慢特征.

SFA算法的优化问题可用如下形式表示:

(1)

\min_{g_{i} (x)}

{\overset{\cdot}{s}}_{i}^{2}

>_t

(2)s.t. <s_i>_t=0

(3)<

s_{i}^{2}

>_t=1

(4)∀i≠j:<s_is_j>_t=0

i,j=1,2,…,m

式中: $s_{i}^{2}$ 为由当前样本计算得到的第i个慢特征统计量; ${\overset{\cdot}{s}}_{i}$ 为慢特征的1阶导数;<·>为时间平均,其定义表达式如下:

(5)<f>_t=

\frac{1}{t_{1} - t_{0}} \int_{t_{0}}^{t_{1}}

f(t)dt

SFA算法通过将正常数据中提取出的慢特征进行线性转化,该过程可表示如下:

(6)s_i=g_i(x)=

{w^{T}}_{i}

i=1,2,…,m

式中:w_i为负载向量.慢特征可表示为

(7)s=Wx

式中:W= ${[w_{1} w_{2} \dots w_{m}]}^{T}$ 为负载矩阵.

在进行SFA算法计算时,首先利用奇异值分解进行白化操作以消除变量之间的相关关系,令R=<x(t)x(t)^T>_t为x(t)的协方差矩阵,R的奇异值分解可写为

(8)R=UΛU^T

式中:Λ为特征对角阵;U为特征矩阵.白化矩阵可以表示为Q=Λ^-1^/²U^Tx,则白化过程可用如下形式表示:

(9)z=Λ^-1^/²U^Tx=Qx

由式(7)和(9)可得:

(10)s=Wx=WQ^-1z=Pz

式中:P=WQ^-1.显然<zz^T>_t=Q<xx^T>Q^T=I,I为单位矩阵且<z>_t=0,又由于约束式(2)和(3)的存在,可以得到如下表达式:

(11)<ss^T>_t=I

则式(11)可写为

(12)<ss^T>_t=P<zz^T>P^T=PP^T=I

由式(12)可知,矩阵P为正交矩阵,则SFA算法的目标是寻找正交矩阵P使得 $< {\overset{\cdot}{s}}_{i}^{2} >_{t}$ 最小化, $< {\overset{\cdot}{s}}_{i}^{2} >_{t}$ 可写为

(13)<

{\overset{\cdot}{s}}_{i}^{2}

>_t=

{p^{T}}_{i}

\overset{\cdot}{z} {\overset{\cdot}{z}}^{T}

>_tp_i

i=1,2,…,m

式中:p_i为特征向量.在实际过程中,样本数据是在离散的时间状态下采集到的,基于时间的导数可以由差分近似计算如下:

(14)

{\overset{\cdot}{z}}_{i}

(t)≈

\frac{z_{i} (t) - z_{i} (t - Δ t)}{Δ t}

i=1,2,…,m

式中:Δt为时间间隔.通过使用协方差矩阵的奇异值分解来解决优化问题,在这种情况下,奇异值分解可以表示为

(15)<

\overset{\cdot}{z} {\overset{\cdot}{z}}^{T}

>_t=P^TΩP

式中:Ω=diag $[λ_{1} λ_{2} \dots λ_{m}]$ 为特征值矩阵,且λ_i升序排列.可以推断出< ${\overset{\cdot}{s}}_{i}^{2}$ >= ${p^{T}}_{i}$ < $\overset{\cdot}{z} {\overset{\cdot}{z}}^{T}$ >_tp_i=λ_i,SFA算法的优化问题是通过选择< ${\overset{\cdot}{s}}_{1}^{2}$ >_t作为最小值,< ${\overset{\cdot}{s}}_{2}^{2}$ >_t作为第二最小值来实现的,以此类推.

SFA负载矩阵的计算可由下式表示:

(16)W=PQ=PΛ^-1^/²U^T

1.2 邻域保持嵌入算法

对于原始数据矩阵X∈Rⁿ^×^m,n为样本个数,NPE通过计算投影矩阵A将X投影到低维空间 Y∈Rⁿ^×^d,其中,d<m,即Y=A^TX.NPE先用k近邻(KNN)准则为每个数据点寻找其k近邻并为其赋权值,定义重构误差为

(17)φ(M)=

\overset{n}{\sum_{i = 1}}

‖x_i-

\sum_{j}

M_ij

x_{i}^{j}

‖²

s.t. $\sum_{j}$ M_ij=1

j=1,2,…,k

式中:x_i为原始样本点;M_ij为样本点x_i与其近邻x_j间的权值;M为权值矩阵;k为近邻个数.

NPE算法的思想是若在高维空间中M_ij可以重构样本x_i,则在低维空间中可通过相同的权值来重构对应低维空间样本点y_i,特征映射可通过最小化目标函数求解:

(18)Φ_NPE=

\sum_{i} {(y_{i} - \sum_{j} M_{ij} y_{i}^{j})}^{2}

s.t. y^Ty= $a_{N}^{T}$ XX^Ta_N=I

i=1,2,…,n;j=1,2,…,k

式中:y_i为原始样本投影到低维空间的样本点,a_N为NPE的投影矩阵A的列向量.a_N可以通过求解下式的广义特征值获得:

(19)XZX^Ta_N=λXX^Ta_N

式中:Z=(I-M)^T(I-M).求解获得最小的d个特征值(λ₁≤λ₂≤…≤λ_d)所对应的特征向量[a_N1a_N2 … a_N_d]^T组成的投影矩阵A.

2 基于GSFA-GNPE动态与静态变量的间歇过程故障检测

2.1 动态变量与静态变量的区分

动态变量是指随时间快速变化的变量,静态变量是指在一段时间内相对恒定的变量.对于原始数据矩阵X= $[x_{1} x_{2} \dots x_{n}]$ ∈Rⁿ^×^m,变量x_i和x_j的互相关系数可以表示为

(20)

ρ_{ij}^{γ}

\frac{1}{n - γ} \frac{\sum_{t}^{=} [x_{i} (t) - μ_{i}] [x_{j} (t + γ) - μ_{j}]}{\sqrt[]{σ_{i}^{2} σ_{j}^{2}}}

i,j=1,2,…,m

式中:μ_i和μ_j为x_i和x_j的平均值; $σ_{i}^{2}$ 和 $σ_{j}^{2}$ 为x_i和x_j的方差;γ为变量间的时间延迟.显然, $ρ_{ij}^{- γ}$ = $ρ_{ji}^{γ}$ ^[19].互相关系数 $ρ_{ij}^{γ}$ 反映了变量i和j的相对方向,一个正的系数表示两个变量以相同的方式发生变化,负相关则表示两个变量以相反的方式发生变化.互相关系数 $ρ_{ij}^{γ}$ 的绝对值反映了变量之间相关性的强弱,绝对值越大,变量间的相关性越强.具有时间延迟γ的x_i顺序相关矩阵 $R_{i}^{γ}$ 可以表示为

(21)

R_{i}^{γ}

[\begin{array}{l} ρ_{i 1}^{- γ} & \dots & ρ_{ij}^{- γ} & \dots & ρ_{im}^{- γ} \\ ρ_{i 1}^{- γ + 1} & \dots & ρ_{ij}^{- λ + 1} & \dots & ρ_{im}^{- γ + 1} \\ ︙ & ︙ & ︙ \\ {\bar{ρ}}_{i 1}^{1} & \dots & {\bar{ρ}}_{ij}^{1} & \dots & {\bar{ρ}}_{im}^{1} \\ ρ_{i 1}^{1} & \dots & ρ_{ij}^{1} & \dots & ρ_{im}^{1} \\ ︙ & ︙ & ︙ \\ ρ_{i 1}^{γ} & \dots & ρ_{ij}^{γ} & \dots & ρ_{im}^{γ} \end{array}]

顺序相关矩阵表示变量和所有变量在不同时滞之间的顺序相关,如果一个变量是静态的,顺序相关矩阵中的所有元素接近于0, $R_{i}^{γ}$ 中的最大绝对值表示为 $ρ_{i}^{\max}$ ,选择一个小的阈值η判断一个变量是否是静态变量,如果 $ρ_{i}^{\max}$ 小于η,则变量i是静态变量,否则变量i是动态变量.

2.2 基于GSFA-GNPE动态与静态变量的间歇过程故障检测

2.2.1 全局慢特征分析算法间歇过程三维数据展开示意图如图1所示,其中:X_i为某一批次的数据矩阵;N为批次数;J为变量数;K为采样时刻数.展开后的数据矩阵X(N×J×K)由一系列的时间切片矩阵X_k(N×J)(k=1,2,…,K)组成,这种展开方式破坏了矩阵X(NK×J)中蕴含的时间序列信息.在GSFA建模时,利用数据近邻点的结构信息构造伪时间序列提供时间序列结构信息.对于样本点x_k_,_i的伪时间序列可表示为

(22)t_k_,_i=

{[\begin{array}{l} x_{k, i} x_{k, i}^{1} & x_{k, i} x_{k, i}^{2} & \dots & x_{k, i} x_{k, i}^{l} \end{array}]}^{T}

式中:x_k_,_i为时间切片矩阵X_k(N×J)中的第i个样本; $x_{k, i}^{c}$ (k=1,2,…,K,c=1,2,…,l)为在对应的时间切片矩阵X_k(N×J)中x_k_,_i的c个近邻.

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 间歇过程三维数据展开示意图

Fig.1 Schematic diagram of 3D data expansion in batch process

更进一步地,可以用如下形式表示时间切片矩阵X_k(N×J)的伪时间序列矩阵:

(23)T_k(2Nl×J)=

{[\begin{array}{l} t_{k, 1} & t_{k, 2} & \dots & t_{k, N} \end{array}]}^{T}

在此基础上,X(NK×J)就被扩展成一个如下式所示的伪时间序列矩阵:

(24)T(2NKl×J)=

[\begin{array}{l} T_{1}^{T} & T_{2}^{T} & \dots & T_{K}^{T} \end{array}]

对于1≤j≤NKl,用τ₂_j_-1和τ₂_j表示伪时间序列T中的第2j-1和第2j个样本点,因此,在样本点τ₂_j_-1和τ₂_j之间可以认为存在时间变化.伪时间序列T的时间变化矩阵ΔT可以表示为

(25)

\begin{array}{l} ΔT = [\begin{array}{l} Δ τ_{1} & Δ τ_{2} & \dots & Δ τ_{NKl} \end{array}] \\ Δ τ_{j} = τ_{2 j} - τ_{2 j - 1} \end{array}\}

SFA的目标函数可用如下形式表示:

(26)Φ_SFA=min

\frac{1}{NKl - 1} {w^{T}}_{j}

\overset{NKl}{\sum_{j = 1}}

(τ₂_j-τ₂_j_-1)×(τ₂_j-τ₂_j_-1)^Tw_j=min

\frac{1}{NKl - 1} {w^{T}}_{j}

ΔTΔT^Tw_j

s.t. $\frac{1}{NK - 1} {w^{T}}_{j} \overset{NK}{\sum_{j = 1}}$ x_j ${x^{T}}_{j}$ w_j=1

式(26)是最小化输出信号的时间变化,由KNN准则构造伪时间序列T(2NKl×J),当在间歇过程处理数据时,SFA 被认为是局部处理方法,因为其仅保留了本地邻域关系而忽略的数据全局结构.为此,提出一种全局慢特征分析算法,在保留邻域信息的基础上分析数据的全局结构.

数据矩阵X全局结构分析的目标是找到一个投影向量ω_j最大化投影样本中的所有样本点y_j= $ω_{j}^{T}$ x_j与其均值之间欧氏距离平方的均值.全局结构目标函数可以表示如下:

(27)Φ_G=max

\frac{1}{NK} \overset{NK}{\sum_{j = 1}}

[y_j-

\bar{y}

]²=max

\frac{1}{NK} \overset{NK}{\sum_{j = 1}}

[

{ω^{T}}_{j}

x_j-

{ω^{T}}_{j} \bar{x}

]²=max

\frac{1}{NK} \overset{NK}{\sum_{j = 1}} {ω^{T}}_{j}

[x_j-

\bar{x}

][x_j-

\bar{x}

]^Tω_j

式中: $\bar{y}$ 为低维空间中的样本均值; $\bar{x}$ = $\frac{1}{NK} \overset{NK}{\sum_{j = 1}}$ x_j为原始数据矩阵中的样本均值.考虑到投影向量ω_j能够写成训练样本的线性组合ω_j=X^Tw_j,且标准化后 $\bar{x}$ =0,则全局目标函数可以写为

(28)

\begin{array}{l} Φ_{G} = \max \frac{1}{NK} \overset{NK}{\sum_{j = 1}} {w^{T}}_{j} x_{j} {x^{T}}_{j} w_{j} = \\ \max {w^{T}}_{j} C w_{j} \\ C = \frac{1}{NK} \overset{NK}{\sum_{j = 1}} x_{j} {x^{T}}_{j} . \end{array}\}

为了保持局部邻域信息和全局结构,GSFA目标函数可以表示为

(29)Φ_GSFA=βΦ_SFA-(1-β)Φ_G=min β

{w^{T}}_{j}

ΔTΔT^Tw_j-max(1-β)

{w^{T}}_{j}

Cw_j=min

{w^{T}}_{j}

[βL-(1-β)C]w_j

j=1,2,…,J

式中:0≤β≤1为松弛因子,用于平衡全局函数和局部函数的比重;L=ΔTΔT^T.考虑到正交约束在建立统计量中的优越性,施加以下正交约束:

${w^{T}}_{j}$ w_j=1

等式中两个目标函数的尺度不同,需要进行标准化,标准化后的目标函数为

(30)min

{w^{T}}_{j}

[

\tilde{β} \tilde{L}

-(1-

\tilde{β}

)

\tilde{C}

]w_j

s.t. ${w^{T}}_{j}$ w_j=1

式中: $\tilde{L}$ =L/ξ(L), $\tilde{C}$ =C/ξ(C)分别为标准化后的L和C;ξ(·)为矩阵的谱半径; $\tilde{β}$ 为标准化后的松弛因子.

求解式(30)的最优化问题等价于求解下式的特征值问题:

(31)[

\tilde{β} \tilde{L}

-(1-

\tilde{β}

)

\tilde{C}

]w_j=

{\tilde{λ}}^{s}_{j}

w_j

式中: ${\tilde{λ}}^{s}_{j}$ 为标准化后第j个慢特征的缓慢变化程度.求解式(31)可得到一系列特征值 ${\tilde{λ}}^{s}_{1}$ ≤ ${\tilde{λ}}^{s}_{2}$ ≤…≤ ${\tilde{λ}}^{s}_{J}$ 及其对应的负载特征向量 $\overset{⌒}{W}$ = $[{\overset{⌒}{w}}_{1} {\overset{⌒}{w}}_{2} \dots {\overset{⌒}{w}}_{J}]$ ,保留前p个广义特征值以构造GSFA负载矩阵 ${\overset{⌒}{W}}_{p}$ = $[{\overset{⌒}{w}}_{1} {\overset{⌒}{w}}_{2} \dots {\overset{⌒}{w}}_{p}]$ .

根据GSFA模型提取前p个慢特征信息后构建S²和平方预测误差(SPE)统计量进行监控,定义如下:

(32)S²=

{s^{T}}_{p}

s_p

(33)SPE=(x-

\hat{x}

)^T(x-

\hat{x}

)

式中: $\hat{x}$ =UΛ^-1^/² ${P^{T}}_{P} {\overset{⌒}{W}}_{P}$ x;P_p和 ${\overset{⌒}{W}}_{P}$ 分别为包含p个慢特征的正交矩阵和负载矩阵.控制限通过核密度估计(KDE)确定^[20].

本文采用的核密度函数表达式如下所示:

(34)ϕ(u)=

\frac{1}{nθ} \overset{NK}{\sum_{j = 1}}

(\frac{u - u_{j}}{θ})

式中:u为采样数据;u_j为数据集中的观测值;θ为核宽度参数; Γ为核函数;选择高斯核函数Γ(j)= $\frac{1}{\sqrt[]{2 π}}$ exp $(- \frac{u^{2}}{2})$ .设置信度α=0.95,S²和SPE的控制限 $S_{α}^{2}$ 和SPE_α计算如下:

(35)

\int_{-}^{\infty}

f(SPE)d(SPE)=

\int_{-}^{\infty} \frac{1}{n} \frac{1}{θ \sqrt[]{2 π}}

\overset{NK}{\sum_{j = 1}}

exp

[- \frac{(SPE - SP E_{j})^{2}}{2 θ^{2}}]

d(SPE)=α

(36)

\int_{-}^{\infty}

f(S²)d(S²)=

\int_{-}^{\infty} \frac{1}{n} \frac{1}{θ \sqrt[]{2 π}}

\overset{NK}{\sum_{j = 1}}

exp

[- \frac{(S^{2} - S_{j}^{2})^{2}}{2 θ^{2}}]

d(S²)=α

式中:SPE_j、S_j为第j个SPE和S²统计量.

2.2.2 GNPE算法 NPE算法在降维过程中只考虑了数据的局部结构,忽略了数据的全局结构信息.为了同时考虑全局和邻域结构,在静态子空间中建立GNPE模型,充分表征故障信息,其中全局结构通过寻找方差最大来保持,其目标函数如下所示:

(37)Φ_G=max

\overset{NK}{\sum_{j}}

‖y_j-

\bar{y}

‖²=max

\overset{NK}{\sum_{j = 1}} a_{G}^{T}

(x_j-

\bar{x}

)(x_j-

\bar{x}

)^Ta_G=max

a_{G}^{T}

Ga_G

j=1,2,…,NK

式中:a_G为全局目标函数的投影向量; $\bar{x}$ = $(\overset{NK}{\sum_{j = 1}} x_{j})$ /(NK);G=(X- $\bar{X}$ )(X- $\bar{X}$ )^T.

GNPE算法的目标函数为

(38)Φ_GNPE=max

\frac{Φ_{G}}{Φ_{NPE}}

= max

\frac{\overset{NK}{\sum_{j}} ‖ y_{j} - \bar{y} ‖^{2}}{\overset{NK}{\sum_{j = 1}} {(y_{j} - \overset{k}{\sum_{i = 1}} M_{ij} y_{i})}^{2}}

引入 ${y^{T}}_{j}$ = ${a^{T}}_{j}$ x_j,式(38)可以转换为

(39)Φ_GNPE=max

\frac{a_{G}^{T} G a_{G}}{a_{N}^{T} R a_{N}}

进一步转换为特征值分解问题:

Gd_i=λ_iRd_i

由最大的g个特征值(λ₁≤λ₂≤…≤λ_g)对应的特征向量[d₁d₂ … d_g]构成投影矩阵D.

GNPE算法将静态空间分为特征空间和残差空间,构造T²和SPE统计量对过程进行监控,其表达式如下所示:

(40)T²=y_new

{\bar{Λ}}_{L}^{- 1}

y_new

(41)SPE=‖x_new-Dy_new‖

(42)y_new=D^TX_new

式中:Λ_L为低维数据矩阵X的协方差矩阵, ${\bar{Λ}}_{L}^{- 1}$ = $\frac{1}{NK - 1} \overset{NK}{\sum_{j = 1}}$ y_j ${y^{T}}_{j}$ ;x_new为测试数据.统计量的控制限根据核密度估计求得.

2.2.3 联合指标建立通过计算相关矩阵,将过程变量分为Ω_D和Ω_S两个子空间,在动态子空间Ω_D进行GSFA建模,计算S²和SPE统计量,在静态子空间Ω_S进行GNPE建模,计算T²和SPE统计量.为了实现对静态子空间和动态子空间的联合监控,通过贝叶斯推断建立一个联合指标,以实现对整个间歇过程的监控.在动态子空间Ω_D中,S²的故障条件概率定义为

(43)

P_{S^{2} (F | Ω_{D})}

\frac{P_{S^{2} (Ω_{D} | F)} P_{S^{2} F}}{P_{S^{2} (Ω_{D} | N)} P_{S^{2} N} + P_{S^{2} (Ω_{D} | F)} P_{S^{2} F}}

(44)

P_{S^{2} (Ω_{D} | N)}

=exp

(- \frac{S_{Ω_{D}}^{2}}{S_{\lim}^{2}})

(45)

P_{S^{2} (X_{Ω_{D}} | F)}

=exp

(- \frac{S_{\lim}^{2}}{S_{Ω_{D}}^{2}})

式中:N和F分别为正常条件和故障条件; $P_{S^{2} (Ω_{D} | N)}$ 和 $P_{S^{2} (X_{Ω_{D}} | F)}$ 为条件概率; $P_{S^{2} N}$ 为正常样本的先验概率, $P_{S^{2} N}$ =α, $P_{S^{2} F}$ =1-α; $S_{\lim}^{2}$ 为S²的控制限; $S_{Ω_{D}}^{2}$ 为动态子空间的S²统计量.

同理,可计算求得 $P_{SPE (F | Ω_{D})}$ 、 $P_{SPE (F | Ω_{S})}$ 和 $P_{T^{2} (F | Ω_{S})}$ ,其中: $P_{SPE (F | Ω_{S})}$ 为动态子空间中SPE的故障概率; $P_{SPE (F | Ω_{S})}$ 为静态子空间中SPE的故障概率; $P_{T^{2} (F | Ω_{S})}$ 为静态子空间中T²的故障概率.动态子空间中S²统计量和静态子空间中T²统计量可建立BIC-C²指标,动态子空间和静态子空间中的SPE统计量可建立BIC-SPE指标,则有:

(46)BIC-C²=

\frac{P_{S^{2} (F | Ω_{D})}^{2} + P_{T^{2} (F | Ω_{S})}^{2}}{P_{S^{2} (F | Ω_{D})} + P_{T^{2} (F | Ω_{S})}}

(47)BIC-SPE=

\frac{P_{SPE (F | Ω_{D})}^{2} + P_{SPE (F | Ω_{S})}^{2}}{P_{SPE (F | Ω_{D})} + P_{SPE (F | Ω_{S})}}

联合指标BIC由如下表达式建立:

(48)BIC=

\frac{P_{SPE (F | Ω_{D})}^{2} + P_{SPE (F | Ω_{S})}^{2} + P_{S^{2} (F | Ω_{D})}^{2} + P_{T^{2} (F | Ω_{S})}^{2}}{P_{SPE (F | Ω_{D})} + P_{SPE (F | Ω_{S})} + P_{S^{2} (F | Ω_{D})} + P_{T^{2} (F | Ω_{S})}}

联合指标BIC-C²、BIC-SPE和BIC的控制限为1-α.

3 基于GSFA-GNPE的动态-静态联合指标间歇过程监控步骤

基于GSFA-GNPE的动态-静态联合指标间歇过程监控方法包括离线建模和在线监测两个阶段,如图2所示.

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 基于GSFA-GNPE的动态-静态联合指标间歇过程监控

Fig.2 Batch process monitoring of dynamic-static joint indicator based on GSFA-GNPE

(1) 离线建模步骤如下.

步骤1 将采集到的正常工况的三维数据矩阵X(N×K×J)展开为X(N×KJ),进行标准化处理后排列成X(NK×J);

步骤2 利用式(21)对每个变量计算顺序相关矩阵 $R_{i}^{γ}$ ;

步骤3 判断 $R_{i}^{γ}$ 中 $ρ_{i}^{\max}$ 是否大于阈值η,若大于η则划分为动态变量,若小于η则划分为静态变量;

步骤4 在划分完成的动态子空间中进行GSFA建模,在静态子空间中进行GNPE建模;

步骤5 计算动态子空间中的控制限 $S_{\lim}^{2}$ 和SPE_lim,计算静态子空间中的控制限 $T_{\lim}^{2}$ 和SPE_lim.

(2) 在线监测步骤如下.

步骤1 对测试样本X_test进行标准化处理;

步骤2 根据离线建模步骤3的结果将过程变量分为动态变量Ω_Dnew和静态变量Ω_Snew;

步骤3 计算动态子空间中的统计量 $S_{new}^{2}$ 和SPE_new,计算静态子空间中的统计量 $T_{new}^{2}$ 和SPE_new;

步骤4 计算联合指标BIC- $C_{new}^{2}$ ,BIC-SPE_new和BIC_new,当联合指标超过其控制限BIC_lim时,表明过程发生故障,否则正常.

4 仿真实验

4.1 数值例子仿真

采用如下数值例子验证所提算法^[21]:

(49)z(i)=

[\begin{array}{l} 0.118 & - 0.191 & 0.27 \\ 0.847 & 0.264 & 0.943 \\ - 0.333 & 0.514 & - 0.217 \end{array}]

z(i-1)+

[\begin{array}{l} 1 & 2 \\ 3 & - 4 \\ - 2 & 1 \end{array}]

u(i-1)

(50)y(i)=z(i)+v(i)

(51)u(i)=

[\begin{array}{l} 0.811 & - 0.226 \\ 0.477 & 0.415 \end{array}]

u(i-1)+

[\begin{array}{l} 0.193 & 0.689 \\ - 0.320 & - 0.749 \end{array}]

h(i-1)

(52)k=

[\begin{array}{l} 0.233 & 0.366 & 0.588 \\ 0.105 & - 0.423 & 0.601 \\ - 0.307 & 0.209 & 0.352 \end{array}]

n(i)+v(i)

式中:h为在 $[- 2,2]$ 上均匀分布的输入向量;输入n= $[n_{1} n_{2} n_{3}]$ 为一个随机向量,n₁、n₂均匀分布在 $[- 1,2]$ 上,n₃均匀分布[0,3]上;v为均值为0、方差为0.1的噪声.输出y、k和u是可测量的,使用8个变量(y₁,y₂,y₃,u₁,u₂,k₁,k₂,k₃)的400个样本进行模型训练.显然,前5个变量是动态的,而后3个变量是静态的,所监测的故障从第200个采样点引入至第400个采样点结束,故障描述如下:

故障1h₂处引入值为2的阶跃故障;

故障2n₁处引入值为3的阶跃故障.

对采集到的400个样本首先对每个变量计算顺序相关矩阵,选择时间延时γ为2^[10].根据统计学理论^[22],一旦两个变量的相关系数绝对值小于0.3,则两个变量不再相关,故本文选择阈值η为0.3.每个顺序相关矩阵中的最大值如表1所示,其中加粗的数据为大于0.3的数值.由表1可以看出,前5个变量顺序相关矩阵中的最大值均大于0.3,被划分为动态变量,其余变量被划分为静态变量,这与所提数值例子变量关系相符.将GSFA-GNPE算法与SFA算法和NPE算法用于数值例子故障检测,置信度 α=0.98.3种算法对故障1和故障2的检测率和误报率如表2所示.由表2可知,GSFA-GNPE算法的检测效果优于SFA算法和NPE算法.

表1 相关系数矩阵的最大值

Tab.1 Max values in cross-correlation matrix

变量	$ρ_{i}^{\max}$	变量	$ρ_{i}^{\max}$
y₁	0.798	u₂	0.798
y₂	0.803	k₁	0.053
y₃	0.600	k₂	0.093
u₁	0.827	k₃	0.100

新窗口打开| 下载CSV

表2 3种算法对故障1和故障2的检测率和误报率

Tab.2 Detection rates and false alarm rates of three algorithms for Fault 1 and Fault 2

故障序号	NPE		SFA		GSFA-GNPE
故障序号	T²	SPE	S²	SPE	BIC-C²	BIC-SPE	BIC
1	0.905/0.190	0.850/0.110	0.795/0.185	0.794/0.160	0.985/0.065	0.920/0.105	0.995/0.095
2	0.75/0.16	0.925/0.085	0.630/0.145	0.63/0.12	0.835/0.145	0.895/0.060	0.945/0.065

新窗口打开| 下载CSV

3种算法对故障1的监控图如图3所示,其中:δ为采样时刻.NPE算法对故障1的T²和SPE监控图如图3(a)所示.从图3(a)中可以看出,NPE算法在采样时刻0~200之间存在大量报警,而0~200时刻为正常时刻,故NPE算法在该阶段存在大量误报警,误报率达19%,且在故障发生时刻200~400之间,有许多采样点统计量在控制限以下,没有及时全面的发出警报.图3(b)为SFA算法对故障1的S²和SPE监控图.从图3(b)中可以看出,在0~200采样点间,SFA算法存在较多的误报,有大量点超过了控制限,且在200~400采样点间有许多漏报.图3(c)为GSFA-GNPE算法的BIC-C²,BIC-SPE和BIC联合指标监控图.从图3(c)中可以看出, BIC-C² 和 BIC-SPE在 0~200 之间只存在少量误报,且在200~400采样点间其检测率均高于SFA算法和NPE算法,BIC联合指标的检测效果最好,不仅其误报率极低,且检测率也达到了99%.综合对比3种算法对于数值例子故障的检测效果,本文所提的GSFA-GNPE算法误报率极小且检测率最大,监控效果最优.

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 3种算法对故障1的监控图

Fig.3 Monitoring chart of Fault 1 by three algorithms

4.2 青霉素发酵过程

本文进一步采用美国伊利诺伊州立理工学院开发的Pensim2.0青霉素发酵过程标准仿真平台得到的过程数据进行仿真验证^[23],该平台可以通过设定不同但都在正常范围内的初始条件,模拟青霉素发酵过程中各变量每个时刻的数据用于分析研究,设定时间为400 h,采样时间为1 h,生产30个批次生长工况数据,并从18个变量中选择10个变量作为监控变量得到数据矩阵X(30×10×400)作为训练样本,如表3所示.故障种类和故障参数设置如表4所示.本文给出5种算法在故障2下的监控图,理想的监控图为在正常阶段(0~200采样点)所有统计量均在控制限以下,在故障阶段(200~400采样点)所有统计量均在控制限以上.

表3 过程变量

Tab.3 Monitoring variables

序号	过程变量	序号	过程变量
1	通风速率	6	溶解氧浓度
2	搅拌功率	7	反应器体积
3	基质馈送率	8	二氧化碳浓度
4	补料温度	9	pH值
5	基质浓度	10	发酵罐温度

新窗口打开| 下载CSV

表4 故障类型

Tab.4 Fault types

故障序号	变量序号	故障类型	幅值/%	引入时间/h
1	1	阶跃	5	200~400
2	2	斜坡	10	200~400
3	2	阶跃	5	200~400
4	1	斜坡	10	200~400

新窗口打开| 下载CSV

在进行GSFA-GNPE建模时,松弛因子β用来平衡数据的局部邻域信息和全局结构信息,作为一个多目标优化问题,β很难找到一个绝对的最优解.因此为了平衡局部信息和全局信息,本文选取β=0.5,相关性阈值η=0.3,时间延迟γ根据文献[24]中的自相关函数法选取γ=17,邻域数k根据经验取k=6,核窗宽控制核函数的径向作用范围,本文通过交叉验证法设置θ=2000.降维维数的确定在过程监控邻域目前尚没有主流的通用方法,本文需将所提算法与DPCA进行比较,因此本文选择与DPCA相同的降维维数,在建立DPCA模型时使用主元贡献度准则确定降维维数,GSFA建模时保留的慢特征数p=6,T²和SPE统计量的统计置信度和贝叶斯推断的置信度均为99%.

根据式(21)计算出顺序相关矩阵,将青霉素发酵过程的10个变量进行动态变量和静态变量的划分,其中变量序号为1、2、4、6、7、8的变量被划分为动态变量,序号为3、5、9、10的变量被划分为静态变量.5种算法对4个故障批次的故障检测率如表5所示.

表5 青霉素发酵过程中4个故障批次故障检测率

Tab.5 Fault detection rates of four failed batches during penicillin fermentation

故障序号	NPE		SFA		TNPE		DPCA		GSFA-GNPE
故障序号	T²	SPE	S²	SPE	T²	SPE	T²	SPE	BIC-C²	BIC-SPE	BIC
1	1/0.085	1/0.005	1/0.05	1/0.005	1/0.06	0.995/0.005	0.99/0.08	1/0.05	1/0.04	1/0.04	1/0.005
2	0.945/0.1	0.95/0	0.875/0.02	0.98/0	0.97/0.15	0.965/0.125	0.915/0.25	0.95/0.135	0.925/0.03	1/0	0.995/0.03
3	0.99/0.08	1/0.01	0.98/0.01	1/0.005	1/0.065	1/0.02	0.99/0.005	1/0.03	1/0.01	1/0	1/0.01
4	0.44/0.105	0.94/0	0.925/0.065	0.475/0	0.65/0.22	0.94/0	0.93/0.07	0.915/0.05	0.815/0.04	0.91/0	0.975/0.04

新窗口打开| 下载CSV

NPE、SFA、TNPE、DPCA和GSFA-GNPE算法在故障2下的监控图如图4所示,故障2为搅拌功率从第200个采样点开始到第400个采样点结束的斜坡故障.由图4(a)可以看出,NPE算法的T²和SPE分别在第210和211个采样点检测出故障,且T²在0~50采样点间存在误报,因此NPE算法对于故障2存在较大的误报及故障检测延时.图4(b)为SFA算法的S²和SPE监控图.从图4(b)中可以看出,S²在第220个采样点才检测出故障发生,且在50~100采样点之间存在较多的误报,SPE在第204个采样点较早检测出故障,且在0~200采样点之间无误报发生.图4(c)为TNPE算法的T²和SPE监控图.从图4(c)中可以看出,TNPE算法的检测延时相较于NPE算法有所降低,原因是TNPE考虑了过程数据的时序特征,但在0~50采样点之间仍存在较多误报.图4(d)为DPCA算法的T²和SPE监控图.从图4(d)中可以看出,DPCA算法检测延时较高且误报也较多,检测效果不佳.图4(e)为GSFA-GNPE算法的BIC-C²,BIC-SPE和BIC联合指标监控图.从图4(e)中可以看出,BIC-C²在第214个采样点检测出故障,且在0~50采样点间有少量故障,而BIC-SPE统计量在故障一发生就立即检测到了故障,且全程无误报漏报,联合指标BIC也在第201个采样点处立即检测出了故障发生,基本不存在检测延迟.通过综合分析5种算法在故障2下的检测率和误报率,GSFA-GNPE算法有着更好的检测效果,原因是GSFA-GNPE算法分别考虑了过程数据的动态特性和静态特性,且针对两种特性分别采取了不同的方法进行监控, 在提取数据局部信息的基础上分析了其全局结构,能够更加全面地捕捉数据特征,因此能够及时有效地反映过程故障.5种算法对4种故障的平均故障检测率F如图5所示.由图5可知,GSFA-GNPE算法的3种指标均有良好的检测效果,进一步说明了本文所提算法的有效性.

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 5种算法对故障2的监控图

Fig.4 Monitoring chart of Fault 2 by five algorithms

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 青霉素发酵过程中4个故障批次的平均故障检测率对比

Fig.5 Comparison of average fault detection rates of four failed batches during penicillin fermentation

5 结语

本文提出一种基于GSFA-GNPE算法的动态-静态联合指标间歇过程监控方法.通过计算顺序相关矩阵,将过程变量划分为动态子空间和静态子空间,并在其中使用对应的算法进行统计分析.在动态子空间中,使用GSFA方法充分提取过程的动态特征,并且通过全局结构分析提取其全局信息特征;在静态子空间中,GNPE方法通过方差最大保持数据的全局结构,同时发挥NPE的局部结构保持能力,能有效提取过程数据的静态信息.然后,使用贝叶斯推断建立联合指标监控整个间歇过程.最后,将本文所提算法应用到数值实例和青霉素发酵仿真过程,研究结果显示所提出的GSFA-GNPE算法具有良好的监控效果.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

WANG

, EDGAR T

, BALDEA

, et al.

A geometric method for batch data visualization, process monitoring and fault detection

[J]. Journal of Process Control, 2018, 67:197-205.

DOI:10.1016/j.jprocont.2017.05.011 URL [本文引用: 1]

[2]

赵春晖, 余万科, 高福荣.

非平稳间歇过程数据解析与状态监控——回顾与展望

[J]. 自动化学报, 2020, 46(10):2072-2091.