CN103853152A - 一种基于ar-pca的间歇过程故障监测方法 - Google Patents

Abstract

本发明提出了一种基于AR-PCA的间歇过程故障监测方法，实现了对具有较强动态性间歇过程的在线监测。传统MPCA方法在监控间歇过程时，未考虑到由于各种随机噪声和干扰的存在而导致变量呈现出相应的自相关和互相关，从而导致在线监控过程中存在大量误报警。本发明首先将测量变量建立多变量自回归(AR)模型，模型系数矩阵采用PLS方法进行辨识，模型阶次采用AIC准则进行辨识，然后对AR模型的残差建立PCA模型，同时本发明算法在线监控新批次数据时引入了训练数据，从而提高了本发明算法的监控效果。本发明弥补了传统MPCA方法在监测具有较强动态性间歇过程时存在大量误报警的不足，对监测实际间歇生产过程具有重大意义。

Description

一种基于AR-PCA的间歇过程故障监测方法

技术领域

本发明是一种间歇过程的故障监测方法，尤其针对具有较强动态性的青霉素发酵过程，应用此方法对生产过程进行故障监测，及时发现生产过程中故障。

背景技术

现代流程工业过程规模不断扩大、复杂性日益增高，投资越来越大，迫使人们对过程生产的安全性和可靠性越来越重视；特别是在一些生物、化工过程中，常常包含高温、高压、易燃、易爆的生产过程，系统一旦发生事故就会造成人员和财产的巨大损失，而且环境污染也要比其他事故严重得多。然而，尽管随着计算机控制技术在工业过程中的逐步普及，许多生产过程实现了自动化，但工业过程中的异常事件管理仍然主要由操作员人工完成。当过程发生异常时，操作员由于受人类自身能力和经验的限制，很容易做出错误判断和行动，这时不但不能使过程恢复正常运行，而且还可能造成更重大的事故。

间歇过程是现代流程工业中常见的一种生产方式，由于其本身具有的灵活性，因此被广泛应用于医药、食品、染料、香料及生化制品等小批量、高附加值产品的生产和制备当中。然而，间歇过程的生产环境及其动态特性与连续过程有明显的区别，频繁地改变生产的产品和工艺操作条件是间歇过程的正常活动方式。间歇过程往往无稳态的工作点，常常从一个稳定状态转为另一个稳定状态，因而可能存在多种状态的组合。间歇过程往往呈现强非线性、动态性和时变特性，其操作复杂度远远大于连续过程，产品质量更容易受到如原材料质量、设备状况、环境条件等不确定性因素的影响。为了提高间歇生产过程与控制系统的可维护性和安全性，并同时提高产品的质量，迫切地需要建立过程监测系统对生产过程进行故障监控与诊断。

目前，以多元统计过程监控为核心的数据驱动方法在流程工业的过程故障监测和诊断领域得到了广泛的关注和研究。然而，传统的多元统计方法在推导相关统计量分布时，假设过程处于稳态，不存在时序相关性，而实际流程工业中的对象几乎很难满足上述条件，当偏离上述假设条件时，采用传统多元统计过程监控算法，就会引起错误的监控结果，导致监控算法的失效。因此，迫切需要一种解决过程动态性的监控方法对生产过程进行监测，进而进行故障诊断。

发明内容

本发明的目的就是针对现阶段过程监测方法的不足，提出了一种间歇过程的故障监测方法。通过引入AR模型对原始测量变量进行预处理，再建立PCA模型进行监测。该方法解决了过程数据具有动态性而导致的监控效果不佳问题，提高了监控性能。

本发明采用了如下的技术方案和实现步骤：

步骤一，使用过程正常运行数据作为数据驱动的训练样本，样本集由同一发酵过程相同工艺下所记录的I批次测量数据构成，X＝(X¹,X²,...,X^I)^T，其中Xⁱ表示第i批次数据；每个批次数据包含K个采样时刻，每个采样时刻采集J个过程变量，即

其中

表示第i批次第k采样时刻采集的数据，

其中表示第i批次中第k采样时刻的第j个过程变量的测量值；

步骤二，将样本集X进行标准化处理；

步骤三，离线建模，分别对标准化处理后的每批次数据建立对应的多变量自回归AR模型，根据建立的AR模型得到每个AR模型的残差，其中第i个AR模型的残差为

L为模型阶次，由赤池信息量准则AIC准则确定，其中，i＝1,2,...,I，

为第i批次数据在k时刻的模型残差，

表示第i批次数据在k时刻的第j个变量的模型残差；

步骤四，将I批次的多边两AR模型的残差E=[E¹,E²,…,E^I]沿变量方向重排列得到E'((K-L)I×J)如下式所示，

$E\′=\left[\begin{array}{cccc}{e}_{L,1}^1& e_{L,2}^1& \·\·\·& e_{L,J}^1\\ e_{L,1}^2& e_{L,2}^2& \·\·\·& e_{L,J}^2\\ \·& & & \\ \·& & & \\ \·& & & \\ e_{L,1}^I& e_{L,2}^I& \·\·\·& e_{L,J}^I\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ e_{K,1}^1& e_{K,2}^1& \·\·\·& e_{K,J}^1\\ e_{K,1}^2& e_{K,2}^2& \·\·\·& e_{K,J}^2\\ \·& & & \\ \·& & & \\ \·& & & \\ e_{K,1}^I& e_{K,2}^I& \·\·\·& e_{K,J}^I\end{array}\right]---\left(1\right)$

E'作为PCA模型的输入，根据主元贡献率>85%确定主元个数R，进而计算出得分矩阵T_PCA与残差矩阵E_PCA，通过得分矩阵与残差矩阵来计算监控统计量Hotelling-T²和SPE以及相应的监控限；

步骤五，在线监测，具体为：在线获取发酵过程k时刻的过程变量，并计算该时刻对应的AR模型残差

公式如下，

$e_k^{\mathrm{new}}=X_k^{\mathrm{new}}-\hat{C}{\hat{X}}_{k-1:k-L}---\left(2\right)$

其中，为I批训练数据AR模型系数矩阵Cⁱ的均值，i＝1,2,…,I;

为在线获取发酵过程k时刻的过程变量；

为训练样本中所有批次I的k时刻之前的L时刻的所有变量的均值构成的向量；

其中，

计算公式如下，

${\hat{X}}_{k-m}=\frac{1}{I-1}\underset{i=1}{\overset{\mathrm{Ii}}{\mathrm{\Σ}}}{X}_{k-m}^i---\left(3\right)$

其中，I为批次数，m=1,2,...,L，是训练样本中第i批次第k-m采样时刻采集的数据；

步骤六，将求得的发酵过程k时刻对应的AR模型的残差

作为PCA模型的输入，得到发酵过程k时刻的Hotelling-T²统计量和SPE统计量，其中所述的PCA模型中的载荷矩阵P_PCA是利用步骤四中所述的PCA模型得到；判断发酵过程k时刻的Hotelling-T²和SPE统计量是否超出步骤四得到的相应控制限，如果二者都未超出，则判定当前过程测量数据正常，重复步骤五、六，直到生产过程结束；否则判定当前过程有故障发生。

步骤二中数据标准化预处理，处理方式如下：

首先计算样本集X的所有时刻上所有过程变量的均值和标准方差，其中第k采样时刻的第j个过程变量的平均值

的计算公式为，表示第i批次中第k采样时刻的第j个过程变量的测量值，k＝1,...,K，j＝1,...,J；第k采样时刻的第j个过程变量的标准方差S_k,j的计算公式为， $S_{k,j}=\sqrt{\frac{1}{I-1}\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{(x_{k,j}^i-{\stackrel{\‾}{X}}_{k,j})}^2},$ k＝1,...,K，j＝1,...,J；

然后对历史数据X进行标准化，其中第i批次中第k采样时刻的第j个过程变量的标准化计算公式如下：

${\tilde{x}}_{k,j}^i=\frac{x_{k,j}^i-{\stackrel{\‾}{X}}_{k,j}}{S_{k,j}}---\left(4\right)$

其中，i＝1,...,I，j＝1,...,J，k＝1,...,K；

步骤三所述的第i批次数据的多变量AR模型如下：

$X_k^i=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{C}_l^i{X}_{k-1}^i+e_k^i---\left(5\right)$

其中，k＝L,L+1,…,K;L通过AIC准则确定，

分别为第i批次数据在k时刻的测量变量和模型残差，为第i批次数据在第k-l时刻的测量变量，

为第i批次数据在第k-l时刻的模型系数向量；第i批次数据的多变量AR模型的系数矩阵

第i批次数据的多变量AR模型的系数矩阵Cⁱ的计算公式如下：

${\left(C^i\right)}^T=W^i{\left({\left(P_{\mathrm{PLS}}^i\right)}^T{W}^i\right)}^{-1}{B}^i{\left(Q_{\mathrm{PLS}}^i\right)}^T---\left(6\right)$

其中，

为第i批次数据的PLS模型输入矩阵的载荷矩阵，第i批次数据的PLS模型输出矩阵的载荷矩阵，Bⁱ为第i批次数据PLS模型的对角回归矩阵，Wⁱ为第i批次数据PLS模型的权重矩阵，T表示转置。

Bⁱ以及Wⁱ是通过对PLS模型进行偏最小二乘辨识得到，其中，PLS模型的输入矩阵为第i批次数据当前时刻k的前L时刻的所有变量 $X_{k-1:k-L}^i\≡[\left(X_{k-1}^i\right),\left(X_{k-2}^i\right),\·\·\·,\left(X_{k-L}^i\right)],$ k＝L,L+1,…,K，输出矩阵为第i批次数据当前时刻k的所有变量

。

有益效果

为了对具有较强动态性的间歇过程进行准确监测，降低误报率及漏报率，本发明提出了一种基于AR-PCA的间歇过程故障监测方法。传统MPCA方法在监控间歇过程时，未考虑到由于各种随机噪声和干扰的存在而导致变量呈现出相应的自相关和互相关，从而导致在线监控过程中存在大量误报警。

首先将测量变量建立多变量自回归(AR)模型，消除了过程动态性的影响，然后对AR模型的残差建立PCA监测模型，对新批次数据进行在线监测。同时本发明算法在线监控新批次数据时引入了训练数据，从而提高了本发明算法的监控效果。

以青霉素发酵过程仿真软件Pensim为实例，进行的实验验证说明，本发明弥补了传统MPCA方法在监测具有较强动态性间歇过程时存在大量误报警的不足，对监测实际间歇生产过程具有重大意义。

附图说明

图1为基于AR-PCA的间歇过程故障监测方法流程图；

图2为基于AR-PCA的间歇过程故障监测方法示意图；

图3为MPCA方法监控正常批次T²统计量结果实验结果，控制限分别为99%和95%；

图4为MPCA方法监控正常批次SPE统计量结果实验结果，控制限分别为99%和95%；

图5为AR-PCA方法监控正常批次T²统计量结果实验结果，控制限分别为99%和95%；

图6为AR-PCA方法监控正常批次SPE统计量结果实验结果，控制限分别为99%和95%；

图7为MPCA方法监控故障批次T²统计量结果实验结果，控制限分别为99%和95%；

图8为MPCA方法监控故障批次SPE统计量结果实验结果，控制限分别为99%和95%；

图9为AR-PCA方法监控故障批次T²统计量结果实验结果，控制限分别为99%和95%；

图10为AR-PCA方法监控故障批次SPE统计量结果实验结果，控制限分别为99%和95%。

具体实施方式

本发明的对象主要是具有较强动态性的间歇过程，应用本发明的方法对生产过程进行故障监测，及时发现生产过程中故障。下面结合青霉素发酵仿真平台Pensim软件，对本发明进一步进行说明。

Pensim仿真平台是由伊利诺科技学院(Illinois Institute ofTechnology，IIT)以Cinar教授为学科带头人的过程建模、监测及控制研究小组于1998－2002年研究开发的。此仿真平台是专门为青霉素发酵过程而设计的，该软件的内核采用基于Bajpai机理模型改进的Birol模型，在此平台上可以简易实现青霉素发酵过程的一系列仿真，相关研究已表明该仿真平台的实用性与有效性，因此已经成为国际上较有影响的青霉素仿真平台。该平台可以对不同操作条件下青霉素生产过程的CO₂浓度、pH值、青霉素浓度及溶氧浓度等进行仿真。通过控制反应过程中的pH值和发酵反应器内的温度，可以使反应在最佳条件下进行。

图1为基于AR-PCA的间歇过程故障监测方法流程图，本发明监测方法是采用离线建模，在线监测的方法。首先基于正常的训练数据离线建立监测模型，然后利用离线建立好的模型在线对新批次数据进行监测；

图2为基于AR-PCA的间歇过程故障监测方法离线建模示意图；

本发明中，离线建模与在线监测的具体实施步骤如下：

步骤一，使用青霉素发酵过程正常运行数据作为数据驱动的训练样本。为了使训练样本数据可靠，同时令训练样本数据足够多，样本集由同一发酵过程相同工艺下所记录的40批次测量数据构成，即I=40，X＝(X¹,X²,...,X⁴⁰)^T，其中Xⁱ表示第i批次数据；每个批次数据包含400个采样时刻，即K=400;每个采样时刻采集10个过程变量进行监测，如表1所示，即J=10，其中

表示第i批次第k采样时刻采集的数据，

其中表示第i批次中第k采样时刻的第j个过程变量的测量值；

表1建立模型所用变量

同时为了使数据更符合实际，每批次的初始条件略有变化，所有测量变量均加入了随机测量噪声。最终建立样本集X(40×(400×10))，可表示为，

$\left[\begin{array}{cccccccccc}{x}_{\mathrm{1,1}}^1& x_{\mathrm{1,2}}^1& \·\·\·& x_{\mathrm{1,10}}^1& \·\·\·& \·\·\·& x_{\mathrm{400,1}}^1& x_{\mathrm{400,2}}^1& \·\·\·& x_{\mathrm{400,10}}^1\\ x_{\mathrm{1,1}}^2& x_{1,2}^2& \·\·\·& x_{\mathrm{1,10}}^2& \·\·\·& \·\·\·& x_{\mathrm{400,1}}^2& x_{\mathrm{400,2}}^2& \·\·\·& x_{\mathrm{400,10}}^2\\ & & \·& & & & & & \·& \\ & & \·& & & & & & \·& \\ & & \·& & & & & & \·& \\ & & \·& & & & & & \·& \\ & & \·& & & & & & \·& \\ & & \·& & & & & & \·& \\ x_{\mathrm{1,1}}^{40}& x_{\mathrm{1,2}}^{40}& \·\·\·& x_{\mathrm{1,10}}^{40}& \·\·\·& \·\·\·& x_{\mathrm{400,1}}^{40}& x_{\mathrm{400,2}}^{40}& \·\·\·& x_{\mathrm{400,10}}^{40}\end{array}\right]---\left(1\right)$

步骤二，对样本集数据进行标准化预处理，包括取数据的中心化、量纲归一化。中心化处理实际是抽取了过程变量在多批次正常运行操作下的平均运行轨迹，这样处理后的数据突显了间歇过程不同操作批次之间的一种正常随机波动，因此可以认为它们近似服从多维正态分布。其数学表达式为，

${\tilde{x}}_{k,j}^i=\frac{x_{k,j}^i-{\stackrel{\‾}{X}}_{k,j}}{S_{k,j}}---\left(2\right)$

其中，其中，i＝1,...,40，j＝1,...,10，k＝1,...,400；

第i批次中第k采样时刻的第j个过程变量的标准化后的数据；为第k采样时刻的第j个过程变量的平均值，计算式为

表示第i批次中第k采样时刻的第j个过程变量的测量值，k＝1,...,400，j＝1,...,10；S_k,j为第k采样时刻的第j个过程变量的标准方差，计算公式为， $S_{k,j}=\sqrt{\frac{1}{I-1}\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{(x_{k,j}^i-{\stackrel{\‾}{X}}_{k,j})}^2}.$

步骤三，离线建模，将标准化处理后的数据，分别对每批次数据建立当前时刻与前L时刻的关系，进而得到批次i的多变量AR模型的残差Eⁱ，每一批次i的多变量AR模型如下，

$X_k^i=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{C}_l^i{X}_{k-1}^i+e_k^i---\left(3\right)$

其中，k＝L,L+1,…,K;L通过AIC准则确定，本发明中经过计算得L=5；

分别为第i批次在k时刻的测量变量和模型残差；

为第i批次在第k-l时刻的测量变量；

为第i批次数据在第k-l时刻的模型系数向量；第i批次数据的多变量AR模型的系数矩阵

考虑到变量的自相关以及互相关的影响，本发明中AR模型的系数矩阵通过PLS的方法进行辨识，PLS模型的输入矩阵为 $X_{k-1:k-L}^i\≡[\left(X_{k-1}^i\right),\left(X_{k-2}^i\right),\·\·\·,\left(X_{k-L}^i\right)],$ 由第i批次某时刻k的前L时刻的所有变量构成；输出矩阵为第i批次某时刻k的所有变量

k＝5,6,…,400,PLS模型定义如下，

$\begin{array}{c}{x}_{k-1:k-L}^i=P_{\mathrm{PLS}}^i{t}_k^i+v_k^i\\ X_l^i=Q_{P\mathrm{LS}}^i{B}^i{t}_k^i+e_k^i\end{array}---\left(4\right)$

其中，

为第i批次前L时刻变量的得分向量；A为主元个数，本发明中A=4；和

分别为第i批次PLS模型中输入矩阵和输出矩阵的载荷矩阵，Bⁱ为第i批次数据PLS模型的对角回归矩阵。

PLS模型建立之后，得分向量按照如下式量可计算：

${\left(t_k^i\right)}^T={\left(X_{k-1:k-L}^i\right)}^T{W}^i{\left({\left(P_{\mathrm{PLS}}^i\right)}^T{W}^i\right)}^{-1}---\left(5\right)$

其中，Wⁱ为第i批次数据PLS模型的权重矩阵，，T表示转置，第i批次AR模型系数矩阵Cⁱ计算公式如下：

${\left(C^i\right)}^T=W^i{\left({\left(P_{\mathrm{PLS}}^i\right)}^T{W}^i\right)}^{-1}{B}^i{\left(Q_{\mathrm{PLS}}^i\right)}^T---\left(6\right)$

本发明中

Bⁱ，Wⁱ均通过非线性迭代最小二乘方法计算得到。

步骤四，用40批次的多变量AR模型的残差构成PCA模型的输入数据，可表示为E＝[E¹,E²,…,E⁴⁰]，将E沿变量方向重排列E'(((400-5)×40)×10)，如下式所示，

$E\′=\left[\begin{array}{cccc}{e}_{5,1}^1& e_{5,2}^1& \·\·\·& e_{5,10}^1\\ e_{5,1}^2& e_{5,2}^2& \·\·\·& e_{5,10}^2\\ \·& & & \\ \·& & & \\ \·& & & \\ e_{51}^{40}& e_{5,2}^{40}& \·\·\·& e_{5,10}^{40}\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ e_{400,1}^1& e_{400,2}^1& \·\·\·& e_{400,10}^1\\ e_{400,1}^2& e_{400,2}^2& \·\·\·& e_{400,10}^2\\ \·& & & \\ \·& & & \\ \·& & & \\ e_{400\£\¬1}^{40}& e_{400,2}^{40}& \·\·\·& e_{400,10}^{40}\end{array}\right]---\left(7\right)$

对重排列后的E'建立PCA模型：

E'＝T_PCA(P_PCA)^T+E_PCA (8)其中，

和

分别为得分矩阵，载荷矩阵，和残差矩阵，根据主元贡献率>85%确定主元个数R=6。

同时，步骤四中，计算监控统计量，其中Hotelling-T²统计量反映的是主元成分趋势变化及幅值偏离模型的程度，SPE统计量反映了测量变量对主元模型偏离度的平方和，然后利用公知常识得到相应的监控限；

步骤五，在线监测，在线获取发酵过程k时刻的过程变量，并计算该时刻对应的AR模型残差为了取得更好的监控效果，对新批次数据建立多变量AR模型时引入训练数据，计算AR模型残差

公式如下，

$e_k^{\mathrm{new}}=X_k^{\mathrm{new}}-\hat{C}{\hat{X}}_{k-1:k-L}---\left(9\right)$

其中，

为40批训练数据AR模型系数矩阵Cⁱ(i＝1,2,…,40)的均值;为新批次k时刻所有变量的测量值；

为训练样本中40批次的某一时刻k对应前L时刻所有变量的均值构成的向量；

其中，

计算公式如下，

${\hat{X}}_{k-m}=\frac{1}{I-1}\underset{i=1}{\overset{\mathrm{Ii}}{\mathrm{\Σ}}}{X}_{k-m}^i---\left(10\right)$

其中，I为批次数，I=40，m=1,2,...，L，

是训练样本中第i批次第k-m采样时刻采集的数据；

在线监控时，计算PCA模型主元和残差公式如下：

$t_{\mathrm{PCA},k}^{\mathrm{new}}=P_{\mathrm{PVA}}^T{\mathrm{\Δe}}_k^{\mathrm{nw}}---\left(11\right)$

$e_{\mathrm{PVA},k}^{\mathrm{nw}}=(I-P_{\mathrm{PCA}}{P}_{\mathrm{PCA}}^T)\mathrm{\Δ}{e}_k^{\mathrm{new}}---\left(12\right)$

其中，

为在线获取发酵过程k时刻AR模型的残差，P_PCA采用步骤四中所述的PCA模型的载荷矩阵P_PCA，Δ为样本集AR模型残差方差的倒数矩阵，由于此时的AR残差已为零均值，此时不需要再次进行中心化。

在线监控时，在线发酵过程k时刻的Htelling-T²定义如下：

$T_k^2=t_{\mathrm{PVA},k}^{\mathrm{new}}{\mathrm{\Λ}}^{-1}{\left(t_{\mathrm{PVA},k}^{\mathrm{new}}\right)}^T---\left(13\right)$

其中，对角阵Λ^-1(R×R)为样本集得分矩阵T_PCA的协方差的逆阵，T²统计量近似服从F分布，当获得了所有时刻的

后，由经验公式得到其监控限计算公式：

$T^2~\frac{R(I^2-1)}{I(I-R)}{F}_{R,I}-R,\mathrm{\α}---\left(14\right)$

其中，R为主元个数R=6；I为批次数I=40；F_R,I-R,α表示自由度为R,I-R的F分布,α为置信度。

同样，在线发酵过程k时刻的SPE统计量定义如下：

${\mathrm{SPE}}_k=e_{\mathrm{PCA},k}^{\mathrm{new}}{\left(e_{\mathrm{PCA}}^{\mathrm{new}}\right)}^T---\left(15\right)$

SPE近似服从χ²分布，其监控控制限计算公式：

$\begin{array}{c}{\mathrm{SPE}}_{\mathrm{k\α}}~g_k{X}_{k,h_k,\mathrm{\α}}^2\\ g_k=\frac{v_k}{2m_k}\\ h_k=\frac{{2m}_k^2}{v_k}\end{array}---\left(16\right)$

其中，式中m_k,v_k分别表示建模数据中所有批次测量数据在k时刻SPE的均值和方差。

为了验证本文研究的监测算法的有效性，与传统MPCA算法进行了比较。分别对同一批次的正常数据和故障数据进行监测。故障批次为由于人为错误操作导致搅拌速率在180小时引幅值为-0.8%阶跃故障直到反应结束。

如图3、4所示，使用传统MPCA方法监控正常批次时，T²虽然没有误报情况，但SPE存在严重的误报，99%与95%的监控限的误报率分别为10%和23%，误报率均超过了1%和5%。图5、6为本发明算法获得的监控结果，T²和SPE均无误报，表明对正常批次监控拥有较好的效果。

图7、8为采用传统MPCA方法对故障批次进行监控的统计量T²和SPE监控结果。可以看出，统计量SPE在前期监控正常数据过程中存在大量误报，并且T²对故障数据的报警时刻严重滞后。图9、10为本发明方法对故障批次进行监控的统计量T²和SPE监控结果。结果表明，T²和SPE对故障数据的报警明显提前的同时，T²和SPE对前期正常数据的误报率均消除。

通过以上两个对比实验，说明本发明的方法是行之有效的。

Claims (5)

1.一种基于AR-PCA的间歇过程故障监测方法，其特征在于,它主要包括下列步骤：

步骤一，使用过程正常运行数据作为数据驱动的训练样本，样本集由同一发酵过程相同工艺下所记录的I批次测量数据构成，X＝(X¹,X²,...,X^I)^T，其中Xⁱ表示第i批次数据；每个批次数据包含K个采样时刻，每个采样时刻采集J个过程变量，即

其中

表示第i批次第k采样时刻采集的数据，

其中

表示第i批次中第k采样时刻的第j个过程变量的测量值；

步骤二，将样本集X进行标准化处理；

步骤三，离线建模，分别对标准化处理后的每批次数据建立对应的多变量自回归AR模型，根据建立的AR模型得到每个AR模型的残差，其中第i个AR模型的残差为

L为模型阶次，由赤池信息量准则AIC准则确定，其中，i＝1,2,...,I，

为第i批次数据在k时刻的模型残差，

表示第i批次数据在k时刻的第j个变量的模型残差；

步骤四，将I批次的多变量AR模型的残差E＝[E¹,E²,…,E^I]沿变量方向重排列得到E'((K-L)I×J)如下式所示，

$E\′=\left[\begin{array}{cccc}{e}_{L,1}^1& e_{L,2}^1& \·\·\·& e_{L,J}^1\\ e_{L,1}^2& e_{L,2}^2& \·\·\·& e_{L,J}^2\\ \·& & & \\ \·& & & \\ \·& & & \\ e_{L,1}^I& e_{L,2}^I& \·\·\·& e_{L,J}^I\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ \·& & & \·\\ e_{K,1}^1& e_{K,2}^1& \·\·\·& e_{K,J}^1\\ e_{K,1}^2& e_{K,2}^2& \·\·\·& e_{K,J}^2\\ \·& & & \\ \·& & & \\ \·& & & \\ e_{K,1}^I& e_{K,2}^I& \·\·\·& e_{K,J}^I\end{array}\right]---\left(1\right)$

E'作为PCA模型的输入，根据主元贡献率>85%确定主元个数R，进而计算出得分矩阵T_PCA与残差矩阵E_PCA，通过得分矩阵与残差矩阵来计算监控统计量Hotelling-T²和SPE以及相应的监控限；

步骤五，在线监测，具体为：在线获取发酵过程k时刻的过程变量，并计算该时刻对应的AR模型残差

公式如下，

$e_k^{\mathrm{new}}=X_k^{\mathrm{new}}-\hat{C}{\hat{X}}_{k-1:k-L}---\left(2\right)$ 其中，

为I批训练数据AR模型系数矩阵Cⁱ的均值，i＝1,2,…,I;

为在线获取发酵过程k时刻的过程变量；

为训练样本中所有批次I的k时刻之前的L时刻的所有变量的均值构成的向量；

其中，

计算公式如下，

${\hat{X}}_{k-m}=\frac{1}{I-1}\underset{i=1}{\overset{\mathrm{Ii}}{\mathrm{\Σ}}}{X}_{k-m}^i---\left(3\right)$

其中，I为批次数，m=1,2,...,L，

是训练样本中第i批次第k-m采样时刻采集的数据；

步骤六，将求得的发酵过程k时刻对应的AR模型的残差

作为PCA模型的输入，得到发酵过程k时刻的Hotelling-T²统计量和SPE统计量，其中所述的PCA模型中的载荷矩阵P_PCA是利用步骤四中所述的PCA模型得到；判断发酵过程k时刻的Hotelling-T²和SPE统计量是否超出步骤四得到的相应控制限，如果二者都未超出，则判定当前过程测量数据正常，重复步骤五、六，直到生产过程结束；否则判定当前过程有故障发生。

2.根据权利要求1所述的一种基于AR-PCA的间歇过程故障监测方法，其特征在于，步骤二中数据标准化预处理，处理方式如下：

首先计算样本集X的所有时刻上所有过程变量的均值和标准方差，其中第k采样时刻的第j个过程变量的平均值

的计算公式为，

表示第i批次中第k采样时刻的第j个过程变量的测量值，k＝1,...,K，j＝1,...,J；第k采样时刻的第j个过程变量的标准方差S_k,j的计算公式为， $S_{k,j}=\sqrt{\frac{1}{I-1}\underset{i=1}{\overset{I}{\mathrm{\Σ}}}{(x_{k,j}^i-{\stackrel{\‾}{X}}_{k,j})}^2},$ k＝1,...,K，j＝1,...,J；

然后对历史数据X进行标准化，其中第i批次中第k采样时刻的第j个过程变量的标准化计算公式如下：

${\tilde{x}}_{k,j}^i=\frac{x_{k,j}^i-{\stackrel{\‾}{X}}_{k,j}}{S_{k,j}}---\left(4\right)$

其中，i＝1,...,I，j＝1,...,J，k＝1,...,K。

3.根据权利要求1所述的一种基于AR-PCA的间歇过程故障监测方法，其特征在于步骤三所述的第i批次数据的多变量AR模型如下：

$X_k^i=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{C}_l^i{X}_{k-1}^i+e_k^i---\left(5\right)$

其中，k＝L,L+1,…,K;L通过AIC准则确定，

分别为第i批次数据在k时刻的测量变量和模型残差，

为第i批次数据在第k-l时刻的测量变量，

为第i批次数据在第k-l时刻的模型系数向量；第i批次数据的多变量AR模型的系数矩阵

4.根据权利要求3所述的一种基于AR-PCA的间歇过程故障监测方法，其特征在于：第i批次数据的多变量AR模型的系数矩阵Cⁱ的计算公式如下：

${\left(C^i\right)}^T=W^i{\left({\left(P_{\mathrm{PLS}}^i\right)}^T{W}^i\right)}^{-1}{B}^i{\left(Q_{\mathrm{PLS}}^i\right)}^T---\left(6\right)$

其中，为第i批次数据的PLS模型输入矩阵的载荷矩阵，第i批次数据的PLS模型输出矩阵的载荷矩阵，Bⁱ为第i批次PLS模型的对角回归矩阵，Wⁱ为第i批次数据PLS模型的权重矩阵，T表示转置。

5.根据权利要求4所述的一种基于AR-PCA的间歇过程故障监测方法，其特征在于：

Bⁱ以及Wⁱ是通过对PLS模型进行偏最小二乘辨识得到，其中，PLS模型的输入矩阵为第i批次数据当前时刻k的前L时刻的所有变量 $X_{k-1:k-L}^i\≡[\left(X_{k-1}^i\right),\left(X_{k-2}^i\right),\·\·\·,\left(X_{k-L}^i\right)],$ k＝L,L+1,…,K，输出矩阵为第i批次数据当前时刻k的所有变量。

Images