CN103970092A - 一种基于自适应fcm的多阶段发酵过程故障监测方法 - Google Patents

Abstract

本发明公开了一种基于自适应FCM的多阶段发酵过程故障监测方法，解决了标准FCM算法用于发酵过程阶段划分时存在的如下问题：无法实现多批次三维数据的聚类；需要人为指定阶段划分个数；随机初始化聚类中心且易受样本噪声及跳变点的影响。具体步骤包括：首先计算每个时刻数据矩阵的相似度指标作为聚类输入样本，依据最大最小聚类法则获取初始聚类中心集合，之后引入聚类有效性函数通过自适应迭代确定最佳聚类个数。本发明实现了基于多个正常操作批次数据的发酵过程阶段划分，使得阶段划分过程更加客观准确，分阶段建模监测模型降低了故障的误报漏报率，对实现发酵过程的控制及故障检测有重要意义。

Description

一种基于自适应FCM的多阶段发酵过程故障监测方法

技术领域

本发明涉及基于MSPM的间歇过程监控与故障诊断技术领域，特别是涉及一种在发酵过程中应用改进的阶段划分算法建立多阶段故障监测模型并实施在线监测的方法。

背景技术

发酵过程是现代流程工业中常见的一种生产方式，被广泛应用于医药、酿酒、生化制品等生产。发酵过程不仅具有一般非线性系统的时变性、大惯性、关联性、不确定性等特点，并且由于发酵过程的中的一些重要参数如菌体浓度和产物浓度等都不可以在线测量，所以发酵过程的控制比一般非线性系统更为复杂。由于发酵过程机理复杂，数据重复性比较差，很难用确定性数学模型来描述，因此，采用统计过程控制等数据驱动技术就成为解决此类复杂非线性生化反应系统的有效途径，已经成为近年来过程控制领域的研究热点之一。发酵过程在线监控和故障监测从历史生产数据出发，通过建立基于数据驱动的统计模型，并用于监视生产过程的进行，及时发现并消除过程的异常状况，实现过程的安全、稳定的运行，最终达到提高产品质量和企业经济效益的目的。

发酵过程是典型的多阶段间歇工业生产过程。目前研究中对发酵过程的在线监控，多采用多向主成分分析(MPCA)和多向偏最小二乘(MPLS)等传统的多元统计方法，这些方法要求变量服从正态分布，利用的仅仅是二阶统计量信息。多向独立成分分析(MICA)是一种基于ICA技术、处理间歇过程三维数据的多元统计方法，不需假设变量满足正态分布，并且利用了信号的高阶统计信息，能更有效地分析处理过程数据，从而能更本质地描述过程特性。

间歇过程的多阶段特性，即间歇过程操作中的过程变量跟随过程操作进程或过程机理特性的变化发生规律性的改变。每个阶段都有其特定的控制目标，有不同的过程主导变量和过程特征，如果将一个完整批次数据作为统计样本建立监控模型，则不能反映出这种特性，监控效果也会很不理想。近几年研究中常用的阶段划分算法即聚类算法，该类算法无需过程先验知识，通过对过程数据进行聚类分析，从而达到分段目的。其中最为典型的有FCM聚类算法，然而该算法用于间歇过程阶段划分时仍存在一些问题尚未解决：无法实现多批次三维数据的聚类；需要事先给定聚类个数，没有明确的度量指标；对初始化的聚类中心十分敏感；易受噪声与孤立点的影响。

发明内容

本发明针对上述现有技术中存在的问题，提出了一种基于自适应FCM的多阶段发酵过程故障监测方法。间歇过程三维数据展开成二维矩阵后，将矩阵间的相似性度量结果作为阶段划分的输入样本，即可实现多批次正常操作三维数据的聚类，又可消除样本噪声与孤立点的影响，增强鲁棒性；依据最大最小距离法则选取初始聚类中心集合，引入聚类有效性函数通过自适应迭代方式确定最佳聚类个数。本发明采用了如下的技术方案和实现步骤：

步骤1：采集发酵过程数据，建立样本集X，所述的样本集由同一发酵过程相同工艺下所记录的I批次测量数据构成，X=(X₁,X₂,...,X_I)^T，其中X_i表示第i批次数据，每个批次包含K个采样时刻，每个采样时刻采集J个过程变量，即X_i=(X_i,1,X_i,2,...,X_i,K)，其中X_i,k表示第i批次第k采样时刻采集的数据，X_i,k=(x_i,k,1,x_i,k,2,...,x_i,k,J)，其中x_i,k,j表示第i批次中第k采样时刻的第j个过程变量的测量值；

步骤2：对样本集X进行标准化处理，处理方式如下：

首先计算样本集X的所有时刻上所有过程变量的均值和标准方差，其中第k采样时刻的第j个过程变量的均值的计算公式为，x_i,k,j表示第i批次中第k采样时刻的第j个过程变量的测量值，k=1,...,K，j=1,...,J；第k采样时刻的第j个过程变量的标准方差s_k,j的计算公式为，k=1,...,K，j=1,...,J；

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

${\tilde{x}}_{i,k,j}=\frac{x_{i,k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}---\left(1\right)$

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

步骤3：将步骤2标准化后的数据重新构造成二维矩阵X'，该矩阵共有J个列向量，即X'=(X₁',X'₂,...,X'_J)，其中第j个列向量X'_j=(X'_j,1,...,X'_j,K)T，X'_j,k=(X'_j,k,1,...,X'_j,k,I)^T，其中X'_j,k,i表示经过步骤2标准化处理后的第j个过程变量第k个采样时刻在第i个批次中对应的值，其中i=1,...,I，j=1,...,J，k=1,...,K；

步骤4：数据矩阵相似性度量

将构造完成的二维矩阵X'中所有I批次上，同一采样时刻采集的所有J个过程变量划分成K个I行J列的时间片矩阵X'_k(I×J),k=1,2,…,K，计算每一个X'_k与其他时刻数据矩阵X_r'的相似度指标D(k，r)，即可得到K×K的相似度矩阵D={D(k，r)k=1，2，…，K；r=1，2，…，K}，作为聚类的输入样本，其中所述的相似度指标D(k，r)计算公式如下：

$D(k,r)=\mathrm{diss}(X_k^{\′},X_r^{\′})=\frac{4}{J}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{({\mathrm{\λ}}_j-0.5)}^2---\left(2\right)$

其中，λ_j为协方差矩阵S_k的第j个特征值，J为过程变量数，

$S_k=\frac{1}{I-1}{P}^T\mathrm{R\; P}---\left(3\right)$

其中，R为矩阵X'_k和X'_r的混合协方差矩阵，P=P_oΛ^-1/2，

$R=\frac{1}{2I-1}(X_{\′}^{\′T}{X}_k^{\′}+X_r^{\′T}{X}_r^{\′})---\left(4\right)$

RP_o=P_oΛ （5）

其中，P_o为对R进行特征值分解得到的正交矩阵，Λ为对角阵，其对角元素为R的特征值；

步骤5：采用最大最小距离法则获取相似度矩阵D的初始聚类中心集合，聚类数为c_D；

步骤6：计算聚类有效性函数的值

令起始聚类数c=2，并在初始聚类中心集合中选取相应个数的聚类中心运行FCM聚类算法，获取当前聚类个数的模糊划分，计算并记录模糊度聚类有效性函数DS的值及其迭代收敛次数；

步骤7：聚类个数c=c+1，重复步骤6，直到c=c_D；自适应选取的最佳聚类个数为DS最小值所对应的聚类个数，最佳聚类个数下得到的模糊划分作为最佳阶段划分；

步骤8：建立离线多阶段监测模型

步骤8.1：根据步骤7得到的最佳阶段划分方式，对步骤3中的重构二维矩阵X'进行阶段划分，获得C个操作子阶段，可表示如下：

$X^{\′}=\left[\begin{array}{c}{X}^{\′\left(1\right)}\\ X^{\′\left(2\right)}\\ \·\\ \·\\ \·\\ X^{\′\left(c\right)}\end{array}\right]=\left[\begin{array}{c}{X}_1^{\′\left(1\right)}\\ \·\\ \·\\ \·\\ X_{k_1}^{\′\left(1\right)}\\ X_{k_1+1}^{\′\left(2\right)}\\ \·\\ \·\\ \·\\ X_{k_2}^{\′\left(2\right)}\\ \·\\ \·\\ \·\\ X_{k_{c-1}}^{\′\left(c\right)}\\ X_{k_{c-1}+1}^{\′\left(c\right)}\\ \·\\ \·\\ \·\\ X_K^{\′\left(c\right)}\end{array}\right]---\left(6\right)$

其中，k₁，k₂，…k_c-1为C个操作子阶段的划分时刻。

步骤8.2：在每个子阶段X'^(t)，t=1,2，…，c内，对过程数据进行PCA白化去相关处理，具体如下：利用主成分分析PCA方法得到白化矩阵Q^t，Q^t=(Λ^t)^-1/2(U^t)^T，其中U^t和Λ^t分别为X'^(t)的协方差矩阵对应的特征向量矩阵和特征值矩阵，之后将X'^(t)进行白化，白化公式为：Z^t=Q^tX'^(t)；

步骤8.3：利用快速ICA算法Fast ICA从Z^t中估计出新的混合矩阵B^t和分离矩阵W^t，再根据S^t=(B^t)^TZ^t得到独立成分S^t，t=1,2，…，c；

步骤8.4：计算每个子阶段的监控统计量(I^t)²和SPE^t，计算公式如下；

(I^t)²=S^t(S^t)^T，t=1，2，…，c （7）

${\mathrm{SPE}}^t=(X^{\′\left(t\right)}-\stackrel{\mathrm{\Λ}}{X^{\′\left(t\right)}}){(X^{\′\left(t\right)}-\stackrel{\mathrm{\Λ}}{X^{\′\left(t\right)}})}^T,t=\mathrm{1,2},\·\·\·,c---\left(8\right)$

其中，S^t为t阶段二维矩阵X'^(t)提取出的独立成分， $\stackrel{\mathrm{\Λ}}{X^{\′\left(t\right)}}=B^t{S}^t,t=\mathrm{1,2},\·\·\·,c$

利用核密度估计法得到每个子阶段监控统计量的置信限；

步骤9：在线监测

步骤9.1：在线获取发酵过程新批次k采样时刻的J个过程变量的数据x_k，并根据步骤2中得到的k时刻的均值和标准方差对其进行标准化得到其中第k采样时刻的第j个过程变量的标准化公式如下：

${\tilde{x}}_{k,j}=\frac{x_{k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}---\left(9\right)$

其中，x_k,j为当前第k采样时刻所采集发酵数据中的第j个过程变量，为第k采样时刻的第j个过程变量的平均值，s_k,j为第k采样时刻的第j个过程变量的标准方差，j=1,...,J，k=1,...,K；

步骤9.2：判断k时刻所在的离线多阶段监测模型，提取标准化后的k时刻采集数据的独立成分s_k，计算公式如下：

$s_k=W^t{\tilde{x}}_k---\left(10\right)$

其中，W^t为多阶段离线建模阶段步骤8.3中所确定的第t个子阶段的分离矩阵；

步骤9.3：计算当前发酵过程k时刻采集数据的监控统计量和SPE_k，计算公式如下：

$I_k^2=s_k{s}_k^T---\left(11\right)$

${\mathrm{SPE}}_k=(\widetilde{x_k}-\stackrel{\stackrel{\mathrm{\Λ}}{~}}{x_k}){(\widetilde{x_k}-\stackrel{\stackrel{\mathrm{\Λ}}{~}}{x_k})}^T---\left(12\right)$

其中，s_k为k时刻采集数据的独立成分，

步骤9.4：将上述计算得到的监控统计量与建模阶段的步骤8.4确定的控制限进行比较，，如果二者都未超出，则判定当前过程测量数据正常，重复步骤9，直到生产过程结束；否则判定当前过程有故障发生。整个算法过程如图1所示。

有益效果

与其他现有技术相比，本发明采用自适应FCM算法实现了发酵过程的自动阶段划分。以数据相似度矩阵作为输入样本，不仅能够有效解决间歇过程三维数据的聚类，亦可消除噪声及离群点对FCM算法的影响；依据最大最小聚类法则获取初始聚类中心，使聚类结果更加准确，减少迭代次数；通过评价函数指导阶段划分个数的最佳选择。本发明方法能客观准确地进行发酵过程的阶段划分，更好地建立多阶段监测模型，从而更有效地进行故障监测，具有一定的实用意义。

附图说明

图1本发明基于自适应FCM的多阶段发酵过程故障监测算法流程图；

图2（a）为聚类有效性函数结果对比图，细实线表示标准FCM算法，虚点线表示本发明改进的FCM算法；

图2（b）为聚类收敛迭代次数对比图，细实线表示标准FCM算法，虚点线表示本发明改进的FCM算法；

图3为阶段划分结果图；

图4（a）为传统MICA方法对阶跃故障批次数据的独立成分平方和I²监测图；

图4（b）为传统MICA方法对阶跃故障批次数据的预测误差平方和SPE监测图；

图5（a）为本发明多阶段MICA方法对阶跃故障批次数据的独立成分平方和I²监测图；

图5（b）为本发明多阶段MICA方法对阶跃故障批次数据的预测误差平方和SPE监测图；

图6（a）为传统MICA方法对斜坡故障批次数据的独立成分平方和I²监测图；

图6（b）为传统MICA方法对斜坡故障批次数据的预测误差平方和SPE监测图；

图7（a）为本发明多阶段MICA方法对斜坡故障批次数据的独立成分平方和I²监测图；

图7（b）为本发明多阶段MICA方法对斜坡故障批次数据的预测误差平方和SPE监测图；

具体实施方式

结合本发明的内容提供如下实施例：

青霉素作为一种抗生素，具有很高的临床医用价值，其生产发酵过程是典型的多阶段间歇操作过程。本文所采用的Pensim仿真平台是由伊利诺科技学院(Illinois Institute of Technology,IIT)的Cinar教授及他所领导的研究技术小组所开发的。此仿真平台是专门为青霉素发酵过程而设计的，具有一定的国际影响力。

影响青霉素发酵过程的可在线观测变量主要有温度、底物流加速度、底物浓度、空气流量、搅拌功率等，本发明选取10个过程变量(表1所示)进行监测。完整的一批青霉素发酵时间大约为400小时，每一小时采样一次，一批可得到400个采样时刻。仿真实验预生成30批次的正常过程数据用于建模，因而可获得三维矩阵X(30×10×400)。

表1建立模型所用变量

对该青霉素发酵过程进行故障监测在计算机上用MATLAB程序实现的具体步骤如下：

步骤1：将上述30批建模数据X(30×10×400)沿批次方向展开得到X(30×(10×400))，可表示为：

$X=\left[\begin{array}{cccccccccc}{x}_{\mathrm{1,1,1}}& x_{\mathrm{1,1,2}}& \·\·\·& x_{\mathrm{1,1,10}}& \·\·\·& \·\·\·& x_{\mathrm{1,400,1}}& x_{\mathrm{1,400,2}}& \·\·\·& x_{\mathrm{1,400,10}}\\ x_{\mathrm{2,1},1}& x_{\mathrm{2,1,2}}& \·\·\·& x_{\mathrm{2,1,10}}& \·\·\·& \·\·\·& x_{\mathrm{2,400,1}}& x_{\mathrm{2,400,2}}& \·\·\·& x_{\mathrm{2,400,10}}\\ & & \·& & & & & & \·& \\ & & \·& & & & & & \·& \\ & & \·& & & & & & \·& \\ & & \·& & & & & & \·& \\ & & \·& & & & & & \·& \\ & & \·& & & & & & \·& \\ x_{\mathrm{30,1},1}& x_{\mathrm{30,1,2}}& \·\·\·& x_{\mathrm{30,1,10}}& \·\·\·& \·\·\·& x_{\mathrm{30,400,1}}& x_{\mathrm{30,400,2}}& \·\·\·& x_{\mathrm{60,400,10}}\end{array}\right]---\left(13\right)$

步骤2：数据标准化处理，包括数据中心化、量纲归一化：

首先按公式计算第k采样时刻的第j个过程变量在所有批次上的均值，其中x_i,k,j为X_30×10×400第i批次中第k采样时刻的第j个过程变量的测量值，k=1,...,400，j=1,...,10；第k采样时刻的第j个过程变量的标准方差s_k,j的计算公式为，k=1,...,400，j=1,...,10；

然后对X_30×10×400进行标准化，其中第i批次中第k采样时刻的第j个过程变量的标准化计算公式如下：

${\tilde{x}}_{k,j}=\frac{x_{k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}---\left(14\right)$

其中，i=1,...,30，j=1,...,10，k=1,...,400；

步骤3：将步骤2标准化后的数据重新构造成二维矩阵X'((30×400)×10)，可表示为：

$X^{\′}=\left[\begin{array}{cccc}{x}_{\mathrm{1,1,1}}& x_{\mathrm{1,1,2}}& \·\·\·& x_{\mathrm{1,1,10}}\\ x_{\mathrm{2,1,1}}& x_{\mathrm{2,1,2}}& \·\·\·& x_{\mathrm{2,1,10}}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ x_{\mathrm{30,1,1}}& x_{\mathrm{30,1,2}}& \·\·\·& x_{\mathrm{30,1,10}}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ x_{\mathrm{1,400,1}}& x_{\mathrm{1.400,2}}& \·\·\·& x_{\mathrm{1,400,10}}\\ x_{\mathrm{2,400,1}}& x_{\mathrm{2,400,2}}& \·\·\·& x_{\mathrm{2,400,10}}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ x_{\mathrm{30,400,1}}& x_{\mathrm{30,400,2}}& \·\·\·& x_{\mathrm{30,400,10}}\end{array}\right]---\left(15\right)$

步骤4：将重构矩阵X'沿时刻方向划分为400个时间片矩阵X'_k(30×10),k=1,2,…,400，计算矩阵相似度得到400×400相似度矩阵D；

步骤5：依据最大最小聚类法则获取相似度矩阵D的初始聚类中心集合 $\mathrm{\Ω}=\left\{\begin{array}{cccc}D(43,:),& D(223,:),& D(2,:),& D(47,:)\\ D(65,:),& D(91,:),& D(315,:)& \end{array}\right\},$ 聚类数c_D=7；

步骤6：令起始聚类个数c=2，在Ω中选取c个聚类中心后，运行FCM算法，得到有效性函数DS的值及相应的迭代次数，如图2所示；

步骤7:聚类个数c=c+1，重复步骤6，直到c=7；自适应选取的最佳聚类个数c=3，最佳聚类个数下得到的模糊划分作为最佳阶段划分，结果如图3所示；

步骤8：分阶段离线建模

步骤8.1：按照步骤7得到的阶段划分方式，对步骤3中的重构二维矩阵X'进行阶段划分，获得三个子阶段：第一阶段为(0-48)小时，第二阶段为(48-213)小时，第三阶段为(213-400)小时；

步骤8.2：在每个子阶段X'^(t),t=1,2,3内，对过程数据进行PCA白化去相关处理，具体如下：利用主成分分析PCA方法得到白化矩阵Q^t，Q^t=(Λ^t)^-1/2(U^t)^T，其中U^t和Λt分别为X'^(t)的协方差矩阵对应的特征向量矩阵和特征值矩阵，之后将X'^(t)进行白化，白化公式为：Z^t=Q^tX'^(t)；

步骤8.3：利用快速ICA算法Fast ICA从Z^t中估计出新的混合矩阵B^t和分离矩阵W^t，再根据S^t=(B^t)^TZ^t得到独立成分S^t,t=1,2,3；

步骤8.4：计算每个子阶段的监控统计量(I^t)²和SPE^t，计算公式如下：

(I^t)²=S^t(S^t)^T，t=1，2，3 （16）

${\mathrm{SPE}}^t=(X^{\′\left(t\right)}-\stackrel{\mathrm{\Λ}}{X^{\′\left(t\right)}}){(X^{\′\left(t\right)}-\stackrel{\mathrm{\Λ}}{X^{\′\left(t\right)}})}^T,t=\mathrm{1,2},\·\·\·,---\left(17\right)$

其中，S^t为t阶段二维矩阵X'^(t)提取出的独立成分，X'^(t)=B^tS^t，t=1，2，3

利用核密度估计法得到每个子阶段监控统计量的置信限；

步骤9：在线监测

步骤9.1：在线获取发酵过程新批次k采样时刻的10个过程变量的数据x_k，并根据步骤2中得到的k时刻的均值和标准方差对其进行标准化得到其中第k采样时刻的第j个过程变量的标准化公式如下：

${\tilde{x}}_{k,j}=\frac{x_{k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}---\left(18\right)$

其中，x_k,j为当前第k采样时刻所采集发酵数据中的第j个过程变量，为第k采样时刻的第j个过程变量的平均值，s_k,j为第k采样时刻的第j个过程变量的标准方差，j=1,...,10，k=1,...,400；

步骤9.2：根据步骤8.1中三个子阶段的时间段范围，判断k采样时刻所处的子阶段，提取标准化后的k时刻采集数据的独立成分s_k，计算公式如下：

$s_k=W^t{\tilde{x}}_k---\left(19\right)$

其中，W^t为多阶段离线建模阶段步骤8.3中所确定的第t个子阶段的分离矩阵；

步骤9.3：计算当前发酵过程k时刻采集数据的监控统计量和SPE_k，计算公式如下：

$I_k^2=s_k{s}_k^T---\left(20\right)$

${\mathrm{SPE}}_k=(\widetilde{x_k}-\stackrel{\stackrel{\mathrm{\Λ}}{~}}{x_k}){(\widetilde{x_k}-\stackrel{\stackrel{\mathrm{\Λ}}{~}}{x_k})}^T---\left(21\right)$

其中，s_k为k时刻采集数据的独立成分，表示采集数据与由独立成分s_k重构回去的估计信号间的残差，

步骤9.4：将上述计算得到的监控统计量与建模阶段的步骤8.4确定的控制限进行比较，，如果二者都未超出，则判定当前过程测量数据正常，重复步骤9，直到生产过程结束；否则判定当前过程有故障发生。

为了验证本文算法的有效性，说明阶段划分结果对监控性能的影响，与传统MICA方法进行了比较。实验选用两批故障批次进行监测：其中一批选为由于人为误操作导致通风速率在200小时引入10%的阶跃故障直至反应结束；另一批选为由于人为因素导致搅拌速率在200小时引入0.2%的斜坡故障直至反应结束。实验结果如图4至图7所示，每幅图分别包括与横坐标平行的线和曲线，其中与横坐标平行的线为通过核密度估计方法确定的控制限，上横线表示99%的控制限，下横线表示95%的控制限，曲线为实时的监测值。如果曲线的值大于控制限的值，说明在此时刻发酵过程发生了故障；否则说明发酵过程运行正常。

图4和图5为分别采用传统MICA方法和本发明提出的多阶段MICA对阶跃故障批次进行监控的效果图。图4、5（a）中与横坐标平行的线为控制限，曲线为实时的I²监测值；图4、5（b）中与横坐标平行的线为控制限，曲线为实时的SPE监测值。两种方法都能及时检测到阶跃故障的发生，但是采用传统MICA方法会有大量误报警。如图4（a）所示的I²统计量在4、9、32、36、38、41时刻处均出现误报警，（b）中SPE统计量的误报率也是高达13.5%；而采用本文提出的多阶段MICA算法图5所示，I²统计量仅有个别误报，99%监控限误报全部消除，SPE统计量的误报率也明显降低。

图6和图7为分别采用传统MICA方法和本文提出方法对斜坡故障批次进行监控的效果图。图6、7（a）中与横坐标平行的线为控制限，曲线为实时的I²监测值；图6、7（b）中与横坐标平行的线为控制限，曲线为实时的SPE监测值。由图6（a）可知，采用传统MICA方法，I²统计量在8、9、10、36、38、41时刻均出现误报警，且在250小时左右才开始报警，滞后时间过长，漏报率高；而图7（a）中本文算法的I²统计量没有出现误报警，214小时处即可对斜坡故障报警。SPE监控统计量也有明显的改善，图6（b）显示，采用传统MICA方法99%与95%监控限的误报率分别为6.5%和3.5%，而在图7（b）中误报率降为1.5%和1%。

通过以上两个实验对比结果，表明本发明的方法能较好地解决间歇过程的阶段划分问题，从而更好地进行故障监测。

Claims (1)

1.一种基于自适应FCM的多阶段发酵过程故障监测方法，其特征在于包括以下步骤：

步骤1：采集发酵过程数据，建立样本集X，所述的样本集由同一发酵过程相同工艺下所记录的I批次测量数据构成，X=(X₁,X₂,...,X_I)^T，其中X_i表示第i批次数据，每个批次包含K个采样时刻，每个采样时刻采集J个过程变量，即X_i=(X_i,1,X_i,2,...,X_i,K)，其中X_i,k表示第i批次第k采样时刻采集的数据，X_i,k=(x_i,k,1,x_i,k,2,...,x_i,k,J)，其中x_i,k,j表示第i批次中第k采样时刻的第j个过程变量的测量值；

步骤2：对样本集X进行标准化处理，处理方式如下：

首先计算样本集X的所有时刻上所有过程变量的均值和标准方差，其中第k采样时刻的第j个过程变量的均值的计算公式为，x_i,k,j表示第i批次中第k采样时刻的第j个过程变量的测量值，k=1,...,K，j=1,...,J；第k采样时刻的第j个过程变量的标准方差s_k,j的计算公式为，k=1,...,K，j=1,...,J；

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

${\tilde{x}}_{i,k,j}=\frac{x_{i,k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}---\left(1\right)$

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

步骤3：将步骤2标准化后的数据重新构造成二维矩阵X'，该矩阵共有J个列向量，即X'=(X₁',X'₂,...,X'_J)，其中第j个列向量X'_j=(X'_j,1,...,X'_j,K)^T，X'_j,k=(X'_j,k,1,...,X'_j,k,I)^T，其中X'_j,k,i表示经过步骤2标准化处理后的第j个过程变量第k个采样时刻在第i个批次中对应的值，其中i=1,...,I，j=1,...,J，k=1,...,K；

步骤4：数据矩阵相似性度量；

将构造完成的二维矩阵X'中所有I批次上，同一采样时刻采集的所有J个过程变量划分成K个I行J列的时间片矩阵X'_k(I×J),k=1,2,…,K，计算每一个X'_k与其他时刻数据矩阵X_r'的相似度指标D(k，r)，即可得到K×K的相似度矩阵D={D(k，r)k=1，2，…，K；r=1，2，…，K}，作为聚类的输入样本，其中所述的相似度指标D(k，r)计算公式如下：

$D(k,r)=\mathrm{diss}(X_k^{\′},X_r^{\′})=\frac{4}{J}\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{({\mathrm{\λ}}_j-0.5)}^2---\left(2\right)$

其中，λ_j为协方差矩阵S_k的第j个特征值，J为过程变量数，

$S_k=\frac{1}{I-1}{P}^T\mathrm{R\; P}---\left(3\right)$

其中，R为矩阵X'_k和X'_r的混合协方差矩阵，P=P_oΛ^-1/2，

$R=\frac{1}{2I-1}(X_{\′}^{\′T}{X}_k^{\′}+X_r^{\′T}{X}_r^{\′})---\left(4\right)$

RP_o=P_oΛ (5)

其中，P_o为对R进行特征值分解得到的正交矩阵，Λ为对角阵，其对角元素为R的特征值；

步骤5：采用最大最小距离法则获取相似度矩阵D的初始聚类中心集合，聚类数为c_D；

步骤6：计算聚类有效性函数的值

令起始聚类数c=2，并在初始聚类中心集合中选取相应个数的聚类中心运行FCM聚类算法，获取当前聚类个数的模糊划分，计算并记录模糊度聚类有效性函数DS的值及其迭代收敛次数；

步骤7：聚类个数c=c+1，重复步骤6，直到c=c_D；自适应选取的最佳聚类个数为DS最小值所对应的聚类个数，最佳聚类个数下得到的模糊划分作为最佳阶段划分；

步骤8：建立离线多阶段监测模型

步骤8.1：根据步骤7得到的最佳阶段划分方式，对步骤3中的重构二维矩阵X'进行阶段划分，获得C个操作子阶段，可表示如下：

$X^{\′}=\left[\begin{array}{c}{X}^{\′\left(1\right)}\\ X^{\′\left(2\right)}\\ \·\\ \·\\ \·\\ X^{\′\left(c\right)}\end{array}\right]=\left[\begin{array}{c}{X}_1^{\′\left(1\right)}\\ \·\\ \·\\ \·\\ X_{k_1}^{\′\left(1\right)}\\ X_{k_1+1}^{\′\left(2\right)}\\ \·\\ \·\\ \·\\ X_{k_2}^{\′\left(2\right)}\\ \·\\ \·\\ \·\\ X_{k_{c-1}}^{\′\left(c\right)}\\ X_{k_{c-1}+1}^{\′\left(c\right)}\\ \·\\ \·\\ \·\\ X_K^{\′\left(c\right)}\end{array}\right]---\left(6\right)$

其中，k₁，k₂，…k_c-1为C个操作子阶段的划分时刻；

步骤8.2：在每个子阶段X'^(t)，t=1,2，…，c内，对过程数据进行PCA白化去相关处理，具体如下：利用主成分分析PCA方法得到白化矩阵Q^t，Q^t=(Λ^t)^-1/2(U^t)^T，其中U^t和Λ^t分别为X'^(t)的协方差矩阵对应的特征向量矩阵和特征值矩阵，之后将X'^(t)进行白化，白化公式为：Z^t=Q^tX'^(t)；

步骤8.3：利用快速ICA算法Fast ICA从Zt中估计出新的混合矩阵B^t和分离矩阵W^t，再根据S^t=(B^t)^TZ^t得到独立成分S^t，t=1,2，…，c；

步骤8.4：计算每个子阶段的监控统计量(I^t)²和SPEt，计算公式如下；

(I^t)²=S^t(S^t)^T，t=1，2，…c (7)

${\mathrm{SPE}}^t=(X^{\′\left(t\right)}-\stackrel{\mathrm{\Λ}}{X^{\′\left(t\right)}}){(X^{\′\left(t\right)}-\stackrel{\mathrm{\Λ}}{X^{\′\left(t\right)}})}^T,t=\mathrm{1,2},\·\·\·,c---\left(8\right)$

其中，S^t为t阶段二维矩阵X'^(t)提取出的独立成分， $\stackrel{\mathrm{\Λ}}{X^{\′\left(t\right)}}=B^t{S}^t,t=\mathrm{1,2},\·\·\·,c$

利用核密度估计法得到每个子阶段监控统计量的置信限；

步骤9：在线监测

步骤9.1：在线获取发酵过程新批次k采样时刻的J个过程变量的数据x_k，并根据步骤2中得到的k时刻的均值和标准方差对其进行标准化得到其中第k采样时刻的第j个过程变量的标准化公式如下：

${\tilde{x}}_{k,j}=\frac{x_{k,j}-{\stackrel{\‾}{x}}_{k,j}}{s_{k,j}}---\left(9\right)$

其中，x_k,j为当前第k采样时刻所采集发酵数据中的第j个过程变量，为第k采样时刻的第j个过程变量的平均值，s_k,j为第k采样时刻的第j个过程变量的标准方差，j=1,...,J，k=1,...,K；

步骤9.2：判断k时刻所在的离线多阶段监测模型，提取标准化后的k时刻采集数据的独立成分s_k，计算公式如下：

$s_k=W^t{\tilde{x}}_k---\left(10\right)$

其中，W^t为多阶段离线建模阶段步骤8.3中所确定的第t个子阶段的分离矩阵；

步骤9.3：计算当前发酵过程k时刻采集数据的监控统计量I_k2和SPE_k，计算公式如下：

$I_k^2=s_k{s}_k^T---\left(11\right)$

${\mathrm{SPE}}_k=(\widetilde{x_k}-\stackrel{\stackrel{\mathrm{\Λ}}{~}}{x_k}){(\widetilde{x_k}-\stackrel{\stackrel{\mathrm{\Λ}}{~}}{x_k})}^T---\left(12\right)$

其中，s_k为k时刻采集数据的独立成分，

步骤9.4：将上述计算得到的监控统计量与建模阶段的步骤8.4确定的控制限进行比较，如果二者都未超出，则判定当前过程测量数据正常，重复步骤9，直到生产过程结束；否则判定当前过程有故障发生。