CN101964021B - 基于递归核主元分析的青霉素发酵过程故障监测方法 - Google Patents

Abstract

一种基于递归核主元分析的青霉素发酵过程故障监测方法，属于故障监测与诊断技术领域，包括：采集通风率，搅拌器功率，基质进给速率，基质进给温度，生成的热量，溶解氧浓度，pH值和二氧化碳浓度；用步骤1标准化后的前N个采样样本，建立初始的监测模型，用RKPCA方法更新模型，并计算特征向量

对连续退火工艺的过程进行故障检测与诊断，当T²统计和SPE统计超出各自的控制限时，认为有故障发生，反之，整个过程正常，本发明方法主要解决数据的非线性和时变性问题，RKPCA通过递归计算训练数据协方差的特征值和特征向量对模型进行更新，结果表明，该方法不仅能大大降低误报警，还提高了故障检测的准确性。

Description

基于递归核主元分析的青霉素发酵过程故障监测方法

技术领域

本发明属于故障监测与诊断技术领域，提出了一种基于递归核主元分析的青霉素发酵过程故障监测方法。

背景技术

间歇及半间歇过程已经广泛应用于化工、发酵、制药、食品生产等诸多领域，青霉素作为抗生素的一种，具有广泛的临床医用价值，而其生产设备是个典型的非线性、动态、多阶段半间歇生产过程。间歇过程操作复杂，操作条件的微小变化、原料中混入杂质等异常状况都会影响到最终产品的产量和质量，如果能对间歇过程在线监测、诊断故障，则能帮助操作人员及时消除故障或暂时停止生产以减小原材料的浪费和设备的占用。

目前，针对间歇过程Nonmikos和Macgregor提出了多元统计监控方法，主要是基于多向主成分分析法(multiway principal component analysis，MPCA)以及多向偏最小二乘法(multiwaypartial least squares，MPLS)，这两种方法都是假设过程变量之间的关系是线性的，它们用于监测简单的间歇过程比较有效，但大多数间歇过程由于操作条件或反应进程的改变，从而具有多阶段性，不同操作阶段的数据动态特性不同，同一操作阶段的变量也往往具有高度非线性，这时用单一的统计模型便不能很好的表现原始数据的信息，还有可能导致一些重要信息的缺失，产生故障的漏报及误报，实际上间歇过程没有稳态操作点，过程变量的轨迹随着时间呈现非线性变化趋势，是典型的非线性动态操作过程，针对间歇过程的非线性特性，Lee等人又提出一种基于核函数的非线性算法：多向核主成分分析(multiway kernel principalcomponent analysis，MKPCA)，提取间歇过程的非线性特征，将间歇过程的三维数据矩阵按纵向展开成二维矩阵并进行标准化，建立过程模型并用于过程的在线监测，核主成分分析(kernel principal component analysis，KPCA)通过非线性核函数把非线性数据映射到高维特征空间，然后在特征空间中执行线性PCA提取特征，KPCA是在高维特征空间中执行PCA，所以无需求解非线性优化问题，并且与其它非线性方法相比，在建模之前它不需要指定主元数目。但是传统KPCA算法存在着不足之处，即KPCA模型是时不变的，而大部分实际工业过程都具有时变的特性，这限制了KPCA在非线性间歇过程中的应用。

发明内容

为弥补上述方法之不足，本发明提出一种基于递归核主元分析的青霉素发酵过程故障监测方法，以达到降低误报警率的目的。

本发明的技术方案是这样实现的：基于递归核主元分析的青霉素发酵过程故障监测方法，包括以下步骤：

步骤1：采集数据，在青霉素发酵过程中，采集的数据包括：通风率，搅拌器功率，基质进给速率，基质进给温度，生成的热量，溶解氧浓度，pH值和二氧化碳浓度数据，其中温度和pH值采用闭环控制，而补料采用开环定值控制，用均值和标准偏差规范采集的数据；

步骤2：计算青霉素发酵过程故障的主要因素P，方法如下：用步骤1标准化后的N个采样样本，建立青霉素发酵过程初始的监测模型，对于青霉素发酵过程中一个新的采样xnew，确定青霉素发酵过程故障的主要因素P和置信限，监测其是否出现异常；如果异常则报警，否则转到步骤3；

其中，青霉素发酵过程故障的主要因素P，公式如下：

$P=\mathrm{\Φ}\left(X\right)\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}B}& \tilde{A}\end{array}\right]U_{\mathrm{\Φ}}^{\′}$

式中，Φ(X)为X＝[x₁，x₂，...，x_N]样本矩阵映射到高维特征空间后的矩阵，N是采样样本个数，青霉素发酵过程初始监测模型主要因素的调节因子

青霉素发酵过程初始模型主要因素的修正矩阵

k(X，x₁)表示X和x₁的内积，K(X)表示样本矩阵的内积，

表示和x₁的内积，

为中间矩阵，

表示中间矩阵的内积，为中间矩阵协方差的特征值阵，U′_Φ为过程变量的特征向量阵；

提取青霉素发酵过程中间矩阵传输因子

表示为：

$\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}B}& \tilde{A}\end{array}\right]=A{\left(U_{\mathrm{\Φ}}^{\′}\right)}^{-1}$

对于青霉素发酵过程，一个新的采样x_new，采集的数据包括：通风率，搅拌器功率，基质进给速率，基质进给温度，生成的热量，溶解氧浓度，pH值和二氧化碳浓度；用初始监测模型得到的主要因素P和置信限，监测新的采样x_new是否出现异常；如果异常则报警，否则转到步骤3；

步骤3：青霉素发酵过程采样x_new是正常数据时，采用递归核主元分析方法更新步骤2的青霉素发酵过程初始的监测模型，并计算更新青霉素发酵过程模型后的青霉素发酵过程故障主要因素公式如下；

$\hat{P}=\mathrm{\Φ}\left(\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]\right)\left[\begin{array}{cc}\tilde{A}& -\frac{1}{h_{\mathrm{\Φ}}^{\′}}\sqrt{\frac{N-1}{N(N-2)}}{B}^{\′}\\ 0^T& \frac{1}{h_{\′\mathrm{\Φ}}^{\′}}\sqrt{\frac{N-1}{N(N-2)}}\end{array}\right]U_{\mathrm{\Φ}}^{\′\′}=\mathrm{\Φ}\left(X_{\mathrm{new}}\right)\hat{A}$

式中，

为更新样本矩阵映射到高维特征空间后的矩阵，青霉素发酵过程更新监测模型主要因素的调节因子

青霉素发酵过程更新模型主要因素的修正矩阵

表示

和x_new的内积；

步骤4：对青霉素发酵过程进行故障监测；

通过Hotelling的T²统计和平方预测误差(SPE)统计进行故障监测，即可判断出青霉素发酵过程是否发生故障，当T²统计和SPE统计超出各自的控制限时，认为有故障发生，反之，整个过程正常，转到步骤3继续更新青霉素发酵过程的监测模型。

步骤2所述的建立青霉素发酵过程初始的监测模型，包括以下步骤：

本发明所提出的RKPCA方法，利用在特征空间中的递归计算来更新样本协方差矩阵的特征值分解。令X＝[x₁，x₂，…，x_N]为青霉素发酵过程样本矩阵，x₁，x₂，…，x_N为青霉素发酵过程采样样本，N是采样样本个数，

为青霉素发酵过程中间矩阵，m是采样变量个数，

为更新后青霉素发酵过程的样本矩阵，x_new为青霉素发酵过程中新样本，它们被映射Φ投影到高维特征空间后分别为Φ(X)，

和Φ(X_new)，因此Φ(X)的均值向量m_Φ和协方差矩阵C^F可计算得到

$m_{\mathrm{\Φ}}=\frac{1}{N}\mathrm{\Φ}\left(\left[\begin{array}{cc}{x}_1& \tilde{X}\end{array}\right]\right)1_N=\frac{1}{N}\mathrm{\Φ}\left(x_1\right)+\frac{N-1}{N}{\tilde{m}}_{\mathrm{\Φ}}---\left(1\right)$

$C^F=\frac{1}{N-1}\stackrel{\‾}{\mathrm{\Φ}}\left(\left[\begin{array}{cc}{x}_1& \tilde{X}\end{array}\right]\right)\stackrel{\‾}{\mathrm{\Φ}}{\left(\left[\begin{array}{cc}{x}_1& \tilde{X}\end{array}\right]\right)}^T$

$=\frac{1}{N-1}{(\mathrm{\Φ}\left(x_1\right)-m_{\mathrm{\Φ}})(\mathrm{\Φ}\left(x_1\right)-m_{\mathrm{\Φ}})}^T+\frac{1}{N-1}\underset{i=2}{\overset{N}{\mathrm{\Σ}}}{(\mathrm{\Φ}\left(x_i\right)-m_{\mathrm{\Φ}})(\mathrm{\Φ}\left(x_i\right)-m_{\mathrm{\Φ}})}^T$

$=\frac{1}{N-1}[\frac{N-1}{N}\mathrm{\Φ}\left(x_1\right)-\frac{N-1}{N}{\tilde{m}}_{\mathrm{\Φ}}]{[\frac{N-1}{N}\mathrm{\Φ}\left(x_1\right)-\frac{N-1}{N}{\tilde{m}}_{\mathrm{\Φ}}]}^T+$

$\frac{1}{N-1}\underset{i=2}{\overset{N}{\mathrm{\Σ}}}[\mathrm{\Φ}\left(x_i\right)-{\tilde{m}}_{\mathrm{\Φ}}+\frac{1}{N}{\tilde{m}}_{\mathrm{\Φ}}-\frac{1}{N}\mathrm{\Φ}\left(x_1\right)]\×{[\mathrm{\Φ}\left(x_i\right)-{\tilde{m}}_{\mathrm{\Φ}}+\frac{1}{N}{\tilde{m}}_{\mathrm{\Φ}}-\frac{1}{N}\mathrm{\Φ}\left(x_1\right)]}^T$

$=\frac{1}{N}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}}){(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})}^T+\frac{1}{N-1}\underset{i=2}{\overset{N}{\mathrm{\Σ}}}(\mathrm{\Φ}\left(x_i\right)-{\tilde{m}}_{\mathrm{\Φ}}){(\mathrm{\Φ}\left(x_i\right)-{\tilde{m}}_{\mathrm{\Φ}})}^T$

$=\frac{1}{N}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}}){(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})}^T+\frac{N-2}{N-1}{\tilde{C}}^F$

$=\frac{N-2}{N-1}\left[\begin{array}{cc}\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& \sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)\end{array}\right]$ (2)

$\×{\left[\begin{array}{cc}\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& \sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)\end{array}\right]}^T$

其中，和

分别表示

的均值向量和协方差矩阵，

为矩阵

投影后去均值的投影矩阵，1_N是由N个1构成的行向量，Φ(x_i)为x_i投影到高维特征空间的值，其中i＝1....N，

为

投影去均值后的矩阵。

令Λ和P分别为C^F的特征值阵和主要因素，

和

分别为

的特征值阵和主要因素。假设

R_Φ是一个正交旋转矩阵，由P＝Φ(X)A，

这里

ξ_i和v_i分别表示

的第i个特征值和特征向量，

ω_i和

分别表示

的第i个特征值和特征向量，将C^F，

分别对角化得到P^TC^FP＝Λ，

由式(2)计算得到其中

$g_{\mathrm{\Φ}}=P^T(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})=A^T[k(X,x_1)-(1/(N-1))K(X,\tilde{X})1_{N-1}].$

令

则

和R_Φ分别对应S_Φ的特征值阵与特征向量阵，由式(2)可得到式(3)

$P^T{C}^{F}P=\frac{1}{N}{P}^T(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}}){(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})}^{T}P+\frac{N-2}{N-1}{P}^T{\tilde{C}}^{F}P$

$=\frac{1}{N}{g}_{\mathrm{\Φ}}{g}_{\mathrm{\Φ}}^T+\frac{N-2}{N-1}{A}^T\mathrm{\Φ}{\left(X\right)}^T\tilde{P}\widetilde{\mathrm{\Λ}}{\tilde{P}}^T\mathrm{\Φ}\left(X\right)A$

$=\frac{1}{N}{g}_{\mathrm{\Φ}}{g}_{\mathrm{\Φ}}^T+\frac{N-2}{N-1}{A}^T\mathrm{\Φ}{\left(X\right)}^T\mathrm{\Φ}\left(\tilde{X}\right)\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^T\mathrm{\Φ}{\left(\tilde{X}\right)}^T\mathrm{\Φ}\left(X\right)A$

$=\frac{1}{N}{g}_{\mathrm{\Φ}}{g}_{\mathrm{\Φ}}^T+\frac{N-2}{N-1}{A}^{T}K(X,\tilde{X})\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^{T}K{(X,\tilde{X})}^{T}A$

$=\mathrm{\Λ}---\left(3\right)$

其中，

表示青霉素发酵过程样本矩阵与中间矩阵的内积；

式(2)中

的奇异值分解满足

$\sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)=\tilde{P}{\widetilde{\mathrm{\Σ}}}_{\mathrm{\Φ}}{\tilde{D}}_{\mathrm{\Φ}}^T---\left(4\right)$

其中，

的主要因素，

为对角阵且满足

为对应的右奇异阵。根据式(4)、式(2)有

$\left[\begin{array}{cc}\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& \sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)\end{array}\right]=$

$\left[\begin{array}{cc}{u}_{\mathrm{\Φ}}& \tilde{P}\end{array}\right]\left[\begin{array}{cc}{h}_{\mathrm{\Φ}}& 0^T\\ \widetilde{\mathrm{\Λ}}{\tilde{P}}^T\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& {\widetilde{\mathrm{\Σ}}}_{\mathrm{\Φ}}\end{array}\right]{\left[\begin{array}{cc}1& 0^T\\ 0_{N-1}& {\tilde{D}}_{\mathrm{\Φ}}\end{array}\right]}^T=---\left(5\right)$

$\left[\begin{array}{cc}{u}_{\mathrm{\Φ}}& \tilde{P}\end{array}\right]\left[\begin{array}{cc}{h}_{\mathrm{\Φ}}& 0^T\\ \widetilde{\mathrm{\Λ}}{R}_{\mathrm{\Φ}}^T{P}^T\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& {\widetilde{\mathrm{\Σ}}}_{\mathrm{\Φ}}\end{array}\right]{\left[\begin{array}{cc}1& 0^T\\ 0_{N-1}& {\tilde{D}}_{\mathrm{\Φ}}\end{array}\right]}^T$

其中，青霉素发酵过程初始监测模型主要因素的调节因子：

$h_{\mathrm{\Φ}}=\left|\right|(I-\tilde{P}\widetilde{\mathrm{\Λ}}{\tilde{P}}^T)\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})\left|\right|$

$=\sqrt{\frac{N-1}{N(N-2)}}\left|\right|(I-\tilde{P}\widetilde{\mathrm{\Λ}}{\tilde{P}}^T)(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})\left|\right|$

$=\sqrt{\frac{N-1}{N(N-2)}}\left|\right|\mathrm{\Φ}\left(x_1\right)-\frac{1}{N-1}\mathrm{\Φ}\left(\tilde{X}\right)1_{N-1}-\mathrm{\Φ}\left(\tilde{X}\right)\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^T(\mathrm{\Φ}{\left(\tilde{X}\right)}^T\mathrm{\Φ}\left(x_1\right)-\frac{1}{N-1}\mathrm{\Φ}{\left(\tilde{X}\right)}^T\mathrm{\Φ}\left(\tilde{X}\right)1_{N-1})\left|\right|$

$=\frac{N-1}{N(N-2)}\left|\right|\mathrm{\Φ}\left(x_1\right)-\frac{1}{N-1}\mathrm{\Φ}\left(\tilde{X}\right)1_{N-1}-\mathrm{\Φ}\left(\tilde{X}\right)\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^T(k\left(\widetilde{X,}{x}_1\right)-\frac{1}{N-1}K\left(\tilde{X}\right)1_{N-1})\left|\right|$

$=\frac{N-1}{N(N-2)}\left|\right|\mathrm{\Φ}\left(x_1\right)-\mathrm{\Φ}\left(\tilde{X}\right)B\left|\right|---\left(6\right)$

$=\frac{N-1}{N(N-2)}\sqrt{1-2B^{T}k(\tilde{X},x_1)+B^{T}K\left(\tilde{X}\right)B}$

$u_{\mathrm{\Φ}}=\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}(I-\tilde{P}\widetilde{\mathrm{\Λ}}{\tilde{P}}^T)(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})$

$=\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}[\mathrm{\Φ}\left(x_1\right)-\mathrm{\Φ}\left(\tilde{X}\right)B]---\left(7\right)$

青霉素发酵过程初始监测模型主要因素的修正矩阵：

$B=\frac{1}{N-1}{1}_{N-1}+\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^T(k(\tilde{X},x_1)-\frac{1}{N-1}K\left(\tilde{X}\right)1_{N-1})---\left(8\right)$

其中

表示青霉素发酵过程中间矩阵的内积，

表示

和x₁的内积；

令

$V_{\mathrm{\Φ}}=\left[\begin{array}{cc}{h}_{\mathrm{\Φ}}& 0^T\\ \widetilde{\mathrm{\Λ}}{R}_{\mathrm{\Φ}}^T{P}^T\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& {\widetilde{\mathrm{\Σ}}}_{\mathrm{\Φ}}\end{array}\right]$

$=\left[\begin{array}{cc}{h}_{\mathrm{\Φ}}& 0^T\\ \widetilde{\mathrm{\Λ}}{R}_{\mathrm{\Φ}}^T{A}^T\sqrt{\frac{N-1}{N(N-2)}}(k(X,x_1)-\frac{1}{N-1}K(X,\tilde{X})1_{N-1})& {\widetilde{\mathrm{\Σ}}}_{\mathrm{\Φ}}\end{array}\right]$

对V_Φ作奇异值分解得到V_Φ＝U′_Φ∑′_ΦD′_ΦT，U′_Φ为特征向量阵，∑′_Φ为对角阵，D′_Φ为对应的右奇异阵，将V_Φ代入式(2)有

$\left[\begin{array}{cc}\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& \sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)\end{array}\right]$

$=\left[\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}\right[\mathrm{\Φ}\left(x_1\right)-\mathrm{\Φ}\left(\tilde{X}\right)B\left]\mathrm{\Φ}\left(\tilde{X}\right)\tilde{A}\right]\×U_{\mathrm{\Φ}}^{\′}{\mathrm{\Σ}}_{\mathrm{\Φ}}^{\′}{D}_{\mathrm{\Φ}}^{\′T}{\left[\begin{array}{cc}1& 0^T\\ 0_{N-1}& {\tilde{D}}_{\mathrm{\Φ}}\end{array}\right]}^T$

$=\mathrm{\Φ}\left(\left[\begin{array}{cc}{x}_1& \tilde{X}\end{array}\right]\right)\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B& \tilde{A}\end{array}\right]\×U_{\mathrm{\Φ}}^{\′}{\mathrm{\Σ}}_{\mathrm{\Φ}}^{\′}{D_{\mathrm{\Φ}}^{\′}}^T{\left[\begin{array}{cc}1& 0^T\\ 0_{N-1}& {\tilde{D}}_{\mathrm{\Φ}}\end{array}\right]}^T---\left(9\right)$

$=\mathrm{\Φ}\left(\left[\begin{array}{cc}\tilde{X}& x_1\end{array}\right]\right)\left[\begin{array}{cc}\tilde{A}& -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B\\ 0^T& \frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}\end{array}\right]\×U_{\mathrm{\Φ}}^{\′}{\mathrm{\Σ}}_{\mathrm{\Φ}}^{\′}{D_{\mathrm{\Φ}}^{\′}}^T{\left[\begin{array}{cc}1& 0^T\\ 0_{N-1}& {\tilde{D}}_{\mathrm{\Φ}}\end{array}\right]}^T$

则有C^F的主要因素P可以表示为

$P=\mathrm{\Φ}\left(\left[\begin{array}{cc}{x}_1& \tilde{X}\end{array}\right]\right)\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B& \tilde{A}\end{array}\right]U_{\mathrm{\Φ}}^{\′}$

$=\mathrm{\Φ}\left(X\right)\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B& \tilde{A}\end{array}\right]U_{\mathrm{\Φ}}^{\′}---\left(10\right)$

又因为P＝Φ(X)A，所以得到式(11)

$A=\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B& \tilde{A}\end{array}\right]U_{\mathrm{\Φ}}^{\′}---\left(11\right)$

根据式(11)，

可由下式计算得到

$\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B& \tilde{A}\end{array}\right]=A{\left(U_{\mathrm{\Φ}}^{\′}\right)}^{-1}$

对于青霉素发酵过程，一个新样本x_new，一个新样本x_new，采集的数据包括：通风率，搅拌器功率，基质进给速率，基质进给温度，生成的热量，溶解氧浓度，pH值和二氧化碳浓度；用步骤2建立初始监测模型得到的主要因素P。青霉素发酵过程在特征空间中的得分向量t∈R^r为

$t=P^T[\mathrm{\Φ}\left(x_{\mathrm{new}}\right)-m_{\mathrm{\Φ}}]=A^T\mathrm{\Φ}{\left(X\right)}^T[\mathrm{\Φ}\left(x_{\mathrm{new}}\right)-\frac{1}{N}\mathrm{\Φ}\left(X\right)1_N]$

$=A^T[k(X,x_{\mathrm{new}})-\frac{1}{N}K\left(X\right)1_N]---\left(12\right)$

其中，P＝[p₁，p₂，…，p_r]，r为保留的核主元数目，k(X，x_new)表示样本矩阵X与新样本x_new的内积向量。新样本x_new的T²和SPE统计量由式(13)和式(14)计算

$T_1^2=t^T{\mathrm{\Λ}}^{-1}t---\left(13\right)$

SPE₁＝[Φ(x_new)-m_Φ]^T(I-PP^T)[Φ(x_new)-m_Φ]                        (14)

其中，Λ为主元的方差阵。T²满足F分布：

$T^2=\frac{r(N^2-1)}{N(N-r)}{F}_{r,N-r}$

其中，N是采样样本的个数，r是保留主元的个数，T²统计量的上限为

$T_{\mathrm{\β}}^2=\frac{r(N^2-1)}{N(N-r)}{F}_{r,N-r,\mathrm{\β}}---\left(15\right)$

其中，β为置信度，而Q统计量的分布满足χ²分布，其控制上限为

Q_β＝gχ²(h)                                        (16)

其中，g＝ρ²/2μ，h＝2μ²/ρ²，μ与ρ²分别对应样本Q统计量的均值和方差，如果

和SPE₁大于它们各自的置信限则报警，说明青霉素发酵过程中发生了异常，否则转到步骤3；

步骤3所述采用递归核主元分析方法更新步骤2的青霉素发酵过程初始的监测模型，并计算更新后青霉素发酵过程模型的主要因素

方法如下：

x_new为一个青霉素发酵过程中新样本并且可用，Φ(x_new)为新样本x_new在特征空间中的投影，

为更新后青霉素发酵过程的样本矩阵，则Φ(X_new)的均值向量和协方差阵

分别为

${\tilde{m}}_{\mathrm{\Φ}}=\frac{1}{N}\mathrm{\Φ}\left(\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]\right)1_N=\frac{N-1}{N}{\tilde{m}}_{\mathrm{\Φ}}+\frac{1}{N}\mathrm{\Φ}\left(x_{\mathrm{new}}\right)---\left(17\right)$

${\hat{C}}^F=\frac{1}{N-1}\stackrel{\‾}{\mathrm{\Φ}}\left(\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]\right)\stackrel{\‾}{\mathrm{\Φ}}{\left(\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]\right)}^T$

$=\frac{N-2}{N-1}\left[\begin{array}{cc}\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_{\mathrm{new}}\right)-{\tilde{m}}_{\mathrm{\Φ}})& \sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)\end{array}\right]\×---\left(18\right)$

${\left[\begin{array}{cc}\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_{\mathrm{new}}\right)-{\tilde{m}}_{\mathrm{\Φ}})& \sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)\end{array}\right]}^T$

由式(2)-(9)我们可以得到：

$V_{\mathrm{\Φ}}^{\′}=\left[\begin{array}{cc}{\widetilde{\mathrm{\Σ}}}_{\mathrm{\Φ}}& \widetilde{\mathrm{\Λ}}{\tilde{A}}^T\sqrt{\frac{N-1}{N(N-2)}}(k(\tilde{X},x_{\mathrm{new}})-\frac{1}{N-1}K\left(\tilde{X}\right)1_{N-1})\\ 0^T& h_{\mathrm{\Φ}}^{\′}\end{array}\right]$

对V′_Φ作奇异值分解得到V′_Φ＝U″_Φ∑″_ΦD″_Φ ^T

进而可以得到

的主要因素

和特征值阵

$\hat{P}=\mathrm{\Φ}\left(\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]\right)\left[\begin{array}{cc}\tilde{A}& -\frac{1}{h_{\mathrm{\Φ}}^{\′}}\sqrt{\frac{N-1}{N(N-2)}}{B}^{\′}\\ 0^T& \frac{1}{h_{\mathrm{\Φ}}^{\′}}\sqrt{\frac{N-1}{N(N-2)}}\end{array}\right]U_{\mathrm{\Φ}}^{\′\′}=\mathrm{\Φ}\left(X_{\mathrm{new}}\right)\hat{A}---\left(19\right)$

$\widehat{\mathrm{\Λ}}=\frac{N-2}{N-1}{{\mathrm{\Σ}}_{\mathrm{\Φ}}^{\′\′}}^2---\left(20\right)$

其中：青霉素发酵过程更新监测模型主要因素的调节因子：

$h_{\mathrm{\Φ}}^{\′}=\frac{N-1}{N(N-2)}\sqrt{1-2B^{\′T}k(\tilde{X},x_{\mathrm{new}})+B^{\′T}K\left(\tilde{X}\right)B^{\′}}---\left(21\right)$

青霉素发酵过程更新监测模型主要因素的修正矩阵：

$B^{\′}=\frac{1}{N-1}{1}_{N-1}+\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^T(k(\tilde{X},x_{\mathrm{new}})-\frac{1}{N-1}K\left(\tilde{X}\right)1_{N-1})---\left(22\right)$

表示青霉素发酵过程中间矩阵和新采样x_new的内积；

步骤4所述的利用Hotelling的T²统计和平方预测误差(SPE)统计进行故障监测，T²统计和平方预测误差(SPE)确定方法如下：

对于青霉素发酵过程中一个新样本z，其在特征空间中的得分向量t∈R^r为

$t={\tilde{P}}^T[\mathrm{\Φ}\left(z\right)-{\tilde{m}}_{\mathrm{\Φ}}]={\hat{A}}^T\mathrm{\Φ}{\left(X_{\mathrm{new}}\right)}^T[\mathrm{\Φ}\left(z\right)-\frac{1}{N}\mathrm{\Φ}\left(X_{\mathrm{new}}\right)1_N]$

$={\hat{A}}^T[k(X_{\mathrm{new}},z)-\frac{1}{N}K\left(X_{\mathrm{new}}\right)1_N]---\left(23\right)$

其中，

r为保留的核主元数目，k(X_new，z)表示青霉素发酵过程更新后的样本矩阵X_new与新样本z的内积向量，青霉素发酵过程新样本z的

和SPE₂统计量由式(24)和式(25)计算

$T_2^2=t^T{\widehat{\mathrm{\Λ}}}^{-1}t---\left(24\right)$

${\mathrm{SPE}}_2={[\mathrm{\Φ}\left(z\right)-{\hat{m}}_{\mathrm{\Φ}}]}^T(I-\hat{P}{\hat{P}}^T)[\mathrm{\Φ}\left(z\right)-{\hat{m}}_{\mathrm{\Φ}}]---\left(25\right)$

其中，为主元的方差阵。

新样本z的和SPE₂统计量的置信限可由式(15)和式(16)得到，若

和SPE₂统计量大于它们的置信限，则认为有故障发生报警，否则转到步骤3；

本发明优点：本发明提出了一种基于递归核主元分析的青霉素发酵过程故障监测方法，主要解决数据的非线性和时变性问题，RKPCA通过递归计算训练数据协方差的特征值和主要因素对模型进行更新，利用此方法进行的过程监测结果表明，该方法不仅能大大降低误报警，还提高了故障监测的准确性。

附图说明

图1为本发明基于递归核主元分析的青霉素发酵过程故障监测方法青霉素发酵过程示意图；

图2为本发明基于递归核主元分析的青霉素发酵过程故障监测方法总流程图；

图3为本发明基于递归核主元分析的青霉素发酵过程故障监测方法基于RKPCA的故障监测流程图；

图4为本发明基于递归核主元分析的青霉素发酵过程监测方法采集到的(a)通风率、(b)搅拌器功率、(c)基质给进速率、(d)基质给进温度、(e)生成的热量、(f)溶解氧的浓度、(g)ph值、(h)二氧化碳浓度示意图；

图5为本发明基于核主元分析的青霉素发酵过程故障监测方法用KPCA得出的青霉素发酵过程的T²统计量；

图6为本发明基于核主元分析的青霉素发酵过程故障监测方法用KPCA得出的青霉素发酵过程的SPE统计量；

图7为本发明基于递归核主元分析的青霉素发酵过程故障监测方法用RKPCA得出的青霉素发酵过程的T²统计量；

图8为本发明基于递归核主元分析的青霉素发酵过程故障监测方法用RKPCA得出的青霉素发酵过程的SPE统计量；

图9为本发明基于递归核主元分析的青霉素发酵过程故障监测方法用RKPCA递归计算的主元个数；

图10为本发明基于递归核主元分析的青霉素发酵过程故障监测方法用RKPCA得出的青霉素发酵过程中斜坡类型故障的T²统计量；

图11为本发明基于递归核主元分析的青霉素发酵过程故障监测方法用RKPCA得出的青霉素发酵过程中斜坡类型故障的SPE统计量。

具体实施方式

下面结合附图和实施例对本发明作进一步详细说明：

青霉素发酵过程是青霉素产生菌在合适的培养基、pH值、温度、空气流量、搅拌等发酵条件下进行生长和合成抗生素的代谢活动。图1为青霉素生产发酵过程示意图，其中，被控变量包括发酵罐的pH值和温度，它们分别通过操纵变量：酸、碱流量和冷、热水流量将其控制在一定的值，主要采用控制器FC控制酸、碱流量和冷、热水的阀门的开度，来调节pH值和温度的。在青霉素发酵过程中温度和pH值采用闭环控制，而补料采用开环定值控制，青霉素发酵的每个批次持续时间是400小时，包含约45小时的预培养阶段和约355小时的间歇补料阶段；

本实施例基于递归核主元分析的青霉素发酵过程故障监测方法，如图2所示，包括以下步骤：

步骤1：采集数据，在青霉素发酵过程中，采集的数据包括：通风率，搅拌器功率，基质进给速率，基质进给温度，生成的热量，溶解氧浓度，pH值和二氧化碳浓度数据，其中温度和pH值采用闭环控制，而补料采用开环定值控制，用均值和标准偏差规范采集的数据；

对于青霉素发酵过程采集两批数据，采用0.5h的采样间隔采集数据，其中一批数据为训练数据用来建立青霉素发酵过程初始的模型，另一批次数据作为测试数据，每批由800个观测数构成，其中在此仿真中取99％的置信限，每个样本包含8个变量，在这里我们给出部分采样数据。分别随机选取了标准数据和实时工况数据中的十组数据如表1和表2所示：

表1.青霉素发酵建模中的十组数据

表2.青霉素发酵测试数据中的十组数据

步骤2：用步骤1标准化后的800个采样样本，建立青霉素发酵过程初始的监测模型，计算青霉素发酵过程故障的主要因素P，并计算置信限；对于青霉素发酵过程中一个新的采样x_new，用初始的监测模型得到的青霉素发酵过程故障主要因素P和置信限，监测其是否出现异常；如果异常则报警，否则转到步骤3；

以青霉素发酵过程800个观测数为样本数据X，以观测数据的后799个数据为中间数据

它们被映射Φ投影到高维特征空间，建立青霉素发酵过程的初始模型，求出中间矩阵的传输因子

根据方程(2)和(10)，计算得到X的协方差矩阵C^F和主要因素P，并求出该模型下的置信限，利用主要因素P可由式(13)和式(14)计算得到青霉素发酵过程中新采样x_new的

和SPE₁的统计量，并它们是否大于各自的置信限，若此时新采样本x_new的统计量没有大于置信限，转到步骤3，本实例通过计算没有发生故障，转到步骤3；

步骤3：在青霉素发酵过程中，青霉素发酵过程新采样x_new是正常数据时，采用递归核主元分析方法更新步骤2的青霉素发酵过程的初始模型，并计算更新后青霉素发酵过程故障的主要因素

x_new为观测数中的一个可以用的新样本，Φ(x_new)为新样本x_new在特征空间中的投影，

为青霉素发酵过程更新后的样本矩阵，则更新后数据样本协方差的

的传输因子和特征值阵

可以有(14)和(15)分别求的，进而可以得到更新后青霉素发酵过程模型后的主要因素

这里分别随机选取了青霉素发酵过程更新模型传输因子

的十组数据如表3所示：

表3.传输因子

十组数据

步骤4：对青霉素发酵过程进行故障监测；

利用上面的两组数据，其中一个批次的数据用来建立KPCA模型，另一批次数据作为测试数据，由图5和图6可见，KPCA由于不具有递归更新的能力，无法跟踪青霉素发酵过程正常工况的变化，进而导致所生成的T²和SPE大约都在450个样本过后超过99％的置信限，产生误报警；

采用递归核主元分析方法更新步骤2的青霉素发酵过程的初始模型，并计算更新青霉素发酵过程模型后的主要因素

为了监控过程，使用T²和SPE统计量，对于一个新样本z∈300，其在特征空间中的得分向量t可有(23)求的，新样本z的T²和SPE的统计量由式(24)和式(25)计算得到。接着根据式(15)和式(16)确定它们的控制限，得到改进的RKPCA对青霉素发酵过程结果如图7和图8所示，从图可以看出，RKPCA由于能够根据新样本来更新模型，可以准确反映青霉素发酵过程正常工况所发生的变化，从而所生成的T²和SPE基本都在99％的置信限之内，在图9中，所保留的主元个数从初始值r＝10变化到最大值r＝13，从而图6和图7中的统计量的置信限也会相应的改变；

为了进一步说明RKPCA对故障的监测能力，青霉素发酵过程的测试数据前700个样本用来展示RKPCA对时变过程的自适应监控效果，后100个样本用作故障数据来进行监测，从第701个样本开始直到第800个样本，基质进给速率以-0.003的斜率下降，从而基质进给速率故障加入到后100个样本中，图9和图10显示了监测结果：RKPCA生成的T²和SPE统计量大约在第770个样本左右监测到加入的故障，此结果表明所提出的RKPCA算法不但能够实现对时变非线性过程的自适应监控，而且能够监测到由故障所引起的过程的重大变化。

Claims (3)

1.基于递归核主元分析的青霉素发酵过程故障监测方法，其特征在于：包括以下步骤：

步骤1：采集数据，在青霉素发酵过程中，采集的数据包括：通风率，搅拌器功率，基质进给速率，基质进给温度，生成的热量，溶解氧浓度，pH值和二氧化碳浓度数据，其中温度和pH值采用闭环控制，而补料采用开环定值控制，用均值和标准偏差标准化采集的数据；

步骤2：计算青霉素发酵过程故障的主要因素P，采用如下方法：用步骤1标准化后的N个采样样本，建立青霉素发酵过程初始的监测模型，对于青霉素发酵过程中一个新的采样x_new，用初始的监测模型得到青霉素发酵过程故障的主要因素P和置信限，监测其是否出现异常；如果异常则报警，否则转到步骤3；

其中，青霉素发酵过程故障的主要因素P，公式如下：

$P=\mathrm{\Φ}\left(X\right)\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B& \tilde{A}\end{array}\right]U_{\mathrm{\Φ}}^{\′}$

式中，Φ(X)为X＝[x₁，x₂，...，x_N]样本矩阵映射到高维特征空间后的矩阵，N是采样样本个数，青霉素发酵过程初始监测模型主要因素的调节因子

$h_{\mathrm{\Φ}}=\frac{N-1}{N(N-2)}\sqrt{1-{2B}^{T}k(X,x_1)+B^{T}k\left(X\right)B},$ 青霉素发酵过程初始模型主要因素的修正矩阵

$B=\frac{1}{N-1}{1}_{N-1}+\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^T(k(\tilde{X},x_1)-\frac{1}{N-1}K\left(\tilde{X}\right)1_{N-1}),$ k(X，x₁)表示X和x₁的内积，K(X)表示样本矩阵的内积，表示

和x₁的内积，

为中间矩阵，

表示中间矩阵的内积，

为中间矩阵协方差的特征值阵，U′_Φ为过程变量的特征向量阵，1_N-1是由N-1个1构成的行向量；

提取青霉素发酵过程中间矩阵传输因子

公式如下：

$\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B& \tilde{A}\end{array}\right]=A{\left(U_{\mathrm{\Φ}}^{\′}\right)}^{-1}$

步骤3：青霉素发酵过程采样X_new是正常数据时，采用递归核主元分析方法更新步骤2的青霉素发酵过程初始的监测模型，并计算更新后的青霉素发酵过程故障的主要因素

公式如下： $\hat{P}=\mathrm{\Φ}\left(\begin{array}{c}\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]\end{array}\right)\left[\begin{array}{cc}\tilde{A}& -\frac{1}{h_{\mathrm{\Φ}}^{\′}}\sqrt{\frac{N-1}{N(N-2)}}{B}^{\′}\\ 0^T& \frac{1}{h_{\mathrm{\Φ}}^{\′}}\sqrt{\frac{N-1}{N(N-2)}}\end{array}\right]U_{\mathrm{\Φ}}^{\′\′}=\mathrm{\Φ}\left(X_{\mathrm{new}}\right)\hat{A}$

式中， $\mathrm{\Φ}\left(X_{\mathrm{new}}\right)=\mathrm{\Φ}\left(\begin{array}{c}\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]\end{array}\right)$ 为更新样本矩阵映射到高维特征空间后的矩阵，青霉素发酵过程更新监测模型主要因素的调节因子 $h_{\mathrm{\Φ}}^{\′}=\frac{N-1}{N(N-2)}\sqrt{1-2B^{\′T}k\left(\tilde{X},x_{\mathrm{new}}\right)+B^{\′T}k\left(\tilde{X}\right)B^{\′}},$ 青霉素发酵过程更新模型主要因素的修正矩阵 $B^{\′}=\frac{1}{N-1}{1}_{N-1}+\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^T(k(\tilde{X},x_{\mathrm{new}})-\frac{1}{N-1}K\left(\tilde{X}\right)1_{N-1}),$

表示

和x_new的内积；

步骤4：对青霉素发酵过程是否发生故障进行监测；

通过监测Hotelling的T²统计和平方预测误差SPE统计进行故障监测，即可判断出青霉素发酵过程是否发生故障，当T²统计和SPE统计超出各自的置信限时，认为有故障发生，反之，整个过程正常，转到步骤3继续更新青霉素发酵过程的监测模型。

2.根据权利要求1所述的基于递归核主元分析的青霉素发酵过程故障监测方法，其特征在于：步骤2所述的建立青霉素发酵过程初始的监测模型，方法如下：

利用在特征空间中的递归计算来更新样本协方差矩阵的特征值分解，令X＝[x₁，x₂，…，x_N]为青霉素发酵过程样本矩阵，x₁，x₂，...，x_N为青霉素发酵过程采样样本，

为青霉素发酵过程中间矩阵，m是采样变量个数， $X_{\mathrm{new}}=\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]$ 为更新后青霉素发酵过程的样本矩阵，x_new为青霉素发酵过程中新样本，它们被映射Φ投影到高维特征空间后分别为Φ(X)，

和Φ(X_new)，因此Φ(X)的均值向量m_Φ和协方差矩阵C^F可计算得到

$m_{\mathrm{\Φ}}=\frac{1}{N}\mathrm{\Φ}\left(\begin{array}{c}\left[\begin{array}{cc}{x}_1& \tilde{X}\end{array}\right]\end{array}\right)1_N=\frac{1}{N}\mathrm{\Φ}\left(x_1\right)+\frac{N-1}{N}{\tilde{m}}_{\mathrm{\Φ}}---\left(1\right)$

$C^F=\frac{1}{N-1}\stackrel{\‾}{\mathrm{\Φ}}\left(\begin{array}{c}\left[\begin{array}{cc}{x}_1& \tilde{X}\end{array}\right]\end{array}\right)\stackrel{\‾}{\mathrm{\Φ}}{\left(\begin{array}{c}\left[\begin{array}{cc}{x}_1& \tilde{X}\end{array}\right]\end{array}\right)}^T$

$=\frac{1}{N-1}(\mathrm{\Φ}\left(x_1\right)-m_{\mathrm{\Φ}}){(\mathrm{\Φ}\left(x_1\right)-m_{\mathrm{\Φ}})}^T+\frac{1}{N-1}\underset{i=2}{\overset{N}{\mathrm{\Σ}}}(\mathrm{\Φ}\left(x_i\right)-m_{\mathrm{\Φ}}){(\mathrm{\Φ}\left(x_i\right)-m_{\mathrm{\Φ}})}^T$

$=\frac{1}{N-1}[\frac{N-1}{N}\mathrm{\Φ}\left(x_1\right)-\frac{N-1}{N}{\tilde{m}}_{\mathrm{\Φ}}]{[\frac{N-1}{N}\mathrm{\Φ}\left(x_1\right)-\frac{N-1}{N}{\tilde{m}}_{\mathrm{\Φ}}]}^T+$

$\frac{1}{N-1}\underset{i=2}{\overset{N}{\mathrm{\Σ}}}[\mathrm{\Φ}\left(x_i\right)-{\tilde{m}}_{\mathrm{\Φ}}+\frac{1}{N}{\tilde{m}}_{\mathrm{\Φ}}-\frac{1}{N}\mathrm{\Φ}\left(x_1\right)]\×{[\mathrm{\Φ}\left(x_i\right)-{\tilde{m}}_{\mathrm{\Φ}}+\frac{1}{N}{\tilde{m}}_{\mathrm{\Φ}}-\frac{1}{N}\mathrm{\Φ}\left(x_1\right)]}^T$

$=\frac{1}{N}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}}){(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})}^T+\frac{1}{N-1}\underset{i=2}{\overset{N}{\mathrm{\Σ}}}(\mathrm{\Φ}\left(x_i\right)-{\tilde{m}}_{\mathrm{\Φ}}){(\mathrm{\Φ}\left(x_i\right)-{\tilde{m}}_{\mathrm{\Φ}})}^T$

$=\frac{1}{N}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}}){(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})}^T+\frac{N-2}{N-1}{\tilde{C}}^F$

$=\frac{N-2}{N-1}\left[\begin{array}{cc}\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& \sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)\end{array}\right]$ (2)

$\×{\left[\begin{array}{cc}\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& \sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)\end{array}\right]}^T$

其中，

和分别表示

的均值向量和协方差矩阵， $\stackrel{\‾}{\mathrm{\Φ}}\left(\begin{array}{c}\left[\begin{array}{cc}{x}_1& \tilde{X}\end{array}\right]\end{array}\right)$ 为矩阵 $\left[\begin{array}{cc}{x}_1& \tilde{X}\end{array}\right]$ 投影后去均值的投影矩阵，1_N是由N个1构成的行向量，Φ(x_i)为x_i投影到高维特征空间的值，其中i＝1...N，

为

投影去均值后的矩阵，

令Λ和P分别为C^F的特征值阵和主要因素，

和分别为

的特征值阵和主要因素，假设

R_Φ是一个正交旋转矩阵，由P＝Φ(X)A，

这里 $A=(I-(1/N)\×E_N)[v_1/\sqrt{{\mathrm{\ξ}}_1}{,v}_2/\sqrt{{\mathrm{\ξ}}_2},\·\·\·,v_i/\sqrt{{\mathrm{\ξ}}_i}],$ ξ_i和v_i分别表示

的第i个特征值和特征向量， $\tilde{A}=(I-(1/(N-1))\×E_{N-1})[{\tilde{v}}_1/\sqrt{{\mathrm{\ω}}_1},{\tilde{v}}_2/\sqrt{{\mathrm{\ω}}_2},\·\·\·,{\tilde{v}}_i/\sqrt{{\mathrm{\ω}}_i}],$ ω_i和

分别表示

的第i个特征值和特征向量，将C^F，

分别对角化得到P^TC^FP＝Λ，

由式(2)计算得到 $[(N-1)/(N-2)]\mathrm{\Λ}-[(N-1)/\left(N(N-2)\right)]g_{\mathrm{\Φ}}{g}_{\mathrm{\Φ}}^T=R_{\mathrm{\Φ}}\widetilde{\mathrm{\Λ}}{R}_{\mathrm{\Φ}}^T,$ 其中 $g_{\mathrm{\Φ}}=P^T(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})=A^T[k(X,x_1)-(1/(N-1))K(x,\tilde{X})1_{N-1}],$

令 $S_{\mathrm{\Φ}}=[(N-1)/(N-2)]\mathrm{\Λ}-[(N-1)/\left(N(N-2)\right)]g_{\mathrm{\Φ}}{g}_{\mathrm{\Φ}}^T,$ 则和R_Φ分别对应S_Φ的特征值阵与特征向量阵，由式(2)可得到式(3)

$P^T{C}^{F}P=\frac{1}{N}{P}^T(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}}){(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})}^{T}P+\frac{N-2}{N-1}{P}^T{\tilde{C}}^{F}P$

$\frac{1}{N}{g}_{\mathrm{\Φ}}{g}_{\mathrm{\Φ}}^T+\frac{N-2}{N-1}{A}^T\mathrm{\Φ}{\left(X\right)}^T\tilde{P}\widetilde{\mathrm{\Λ}}{\tilde{P}}^T\mathrm{\Φ}\left(X\right)A$

$=\frac{1}{N}{g}_{\mathrm{\Φ}}{g}_{\mathrm{\Φ}}^T+\frac{N-2}{N-1}{A}^T\mathrm{\Φ}{\left(X\right)}^T\mathrm{\Φ}\left(\tilde{X}\right)\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^T\mathrm{\Φ}{\left(\tilde{X}\right)}^T\mathrm{\Φ}\left(X\right)A$

$=\frac{1}{N}{g}_{\mathrm{\Φ}}{g}_{\mathrm{\Φ}}^T+\frac{N-2}{N-1}{A}^{T}K(X,\tilde{X})\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^{T}K{(X,\tilde{X})}^{T}A---\left(3\right)$

$=\mathrm{\Λ}$

其中，

表示青霉素发酵过程样本矩阵与中间矩阵的内积；

式(2)中

的奇异值分解满足

$\sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)=\tilde{P}{\widetilde{\mathrm{\Σ}}}_{\mathrm{\Φ}}{\tilde{D}}_{\mathrm{\Φ}}^T---\left(4\right)$

其中，

为

的主要因素，

为对角阵且满足

为对应的右奇异阵，根据式(4)、式(2)有

$\left[\begin{array}{cc}\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& \sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)\end{array}\right]=$

$\left[\begin{array}{cc}{u}_{\mathrm{\Φ}}& \tilde{P}\end{array}\right]\left[\begin{array}{cc}{h}_{\mathrm{\Φ}}& 0^T\\ \widetilde{\mathrm{\Λ}}{\tilde{P}}^T\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& {\widetilde{\mathrm{\Σ}}}_{\mathrm{\Φ}}\end{array}\right]{\left[\begin{array}{cc}1& 0^T\\ 0_{N-1}& {\tilde{D}}_{\mathrm{\Φ}}\end{array}\right]}^T=---\left(5\right)$

$\left[\begin{array}{cc}{u}_{\mathrm{\Φ}}& \tilde{P}\end{array}\right]\left[\begin{array}{cc}{h}_{\mathrm{\Φ}}& 0^T\\ \widetilde{\mathrm{\Λ}}{R}_{\mathrm{\Φ}}^T{P}^T\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& {\widetilde{\mathrm{\Σ}}}_{\mathrm{\Φ}}\end{array}\right]{\left[\begin{array}{cc}1& 0^T\\ 0_{N-1}& {\tilde{D}}_{\mathrm{\Φ}}\end{array}\right]}^T$

其中，青霉素发酵过程初始监测模型主要因素的调节因子：

$h_{\mathrm{\Φ}}=\left|\right|(I-\tilde{P}\widetilde{\mathrm{\Λ}}{\tilde{P}}^T)\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})\left|\right|$

$=\sqrt{\frac{N-1}{N(N-2)}}\left|\right|(I-\tilde{P}\widetilde{\mathrm{\Λ}}{\tilde{P}}^T)(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})\left|\right|$

$=\sqrt{\frac{N-1}{N(N-2)}}\left|\right|\mathrm{\Φ}\left(x_1\right)-\frac{1}{N-1}\mathrm{\Φ}\left(\tilde{X}\right)1_{N-1}-\mathrm{\Φ}\left(\tilde{X}\right)\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^T(\mathrm{\Φ}{\left(\tilde{X}\right)}^T\mathrm{\Φ}\left(x_1\right)-\frac{1}{N-1}\mathrm{\Φ}{\left(\tilde{X}\right)}^T\mathrm{\Φ}\left(\tilde{X}\right)1_{N-1})\left|\right|$

$=\frac{N-1}{N(N-2)}\left|\right|\mathrm{\Φ}\left(x_1\right)-\frac{1}{N-1}\mathrm{\Φ}\left(\tilde{X}\right)1_{N-1}-\mathrm{\Φ}\left(\tilde{X}\right)\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^T(k(\tilde{X},x_1)-\frac{1}{N-1}K\left(\tilde{X}\right)1_{N-1})\left|\right|$

$=\frac{N-1}{N(N-2)}\left|\right|\mathrm{\Φ}\left(x_1\right)-\mathrm{\Φ}\left(\tilde{X}\right)B\left|\right|$

$=\frac{N-1}{N(N-2)}\sqrt{1-2B^{T}k(\tilde{X},x_1)+B^{T}k\left(\tilde{X}\right)B}$ (6)

$u_{\mathrm{\Φ}}=\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}(I-\tilde{P}\widetilde{\mathrm{\Λ}}{\tilde{P}}^T)(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})$

$=\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}[\mathrm{\Φ}\left(x_1\right)-\mathrm{\Φ}\left(\tilde{X}\right)B]---\left(7\right)$

青霉素发酵过程初始监测模型主要因素的修正矩阵：

$B=\frac{1}{N-1}{1}_{N-1}+\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^T(k(\tilde{X},x_1)-\frac{1}{N-1}K\left(\tilde{X}\right)1_{N-1})---\left(8\right)$

其中

表示青霉素发酵过程中间矩阵的内积，

表示

和x₁的内积；

令 $V_{\mathrm{\Φ}}=\left[\begin{array}{cc}{h}_{\mathrm{\Φ}}& 0^T\\ \widetilde{\mathrm{\Λ}}{R}_{\mathrm{\Φ}}^T{P}^T\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& {\widetilde{\mathrm{\Σ}}}_{\mathrm{\Φ}}\end{array}\right]$

$=\left[\begin{array}{cc}{h}_{\mathrm{\Φ}}& 0^T\\ \widetilde{\mathrm{\Λ}}{R}_{\mathrm{\Φ}}^T{A}^T\sqrt{\frac{N-1}{N(N-2)}}(k(X,x_1)-\frac{1}{N-1}K(X,\tilde{X})1_{N-1})& {\widetilde{\mathrm{\Σ}}}_{\mathrm{\Φ}}\end{array}\right]$

对V_Φ作奇异值分解得到V_Φ＝U_Φ∑＇_ΦD′_Φ ^T，U′_Φ为特征向量阵，∑′_Φ为对角阵，D′_Φ为对应的右奇异阵，将V_Φ代入式(2)有

$\left[\begin{array}{cc}\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_1\right)-{\tilde{m}}_{\mathrm{\Φ}})& \sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)\end{array}\right]$

$=\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}[\mathrm{\Φ}\left(x_1\right)-\mathrm{\Φ}\left(\tilde{X}\right)B]& \mathrm{\Φ}\left(\tilde{X}\right)\tilde{A}\end{array}\right]\×U_{\mathrm{\Φ}}^{\′}{\mathrm{\Σ}}_{\mathrm{\Φ}}^{\′}{D_{\mathrm{\Φ}}^{\′}}^T{\left[\begin{array}{cc}1& 0^T\\ 0_{N-1}& {\tilde{D}}_{\mathrm{\Φ}}\end{array}\right]}^T$

$=\mathrm{\Φ}\left(\begin{array}{c}\left[\begin{array}{cc}{x}_1& \tilde{X}\end{array}\right]\end{array}\right)\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B& \tilde{A}\end{array}\right]\×U_{\mathrm{\Φ}}^{\′}{\mathrm{\Σ}}_{\mathrm{\Φ}}^{\′}{D_{\mathrm{\Φ}}^{\′}}^T{\left[\begin{array}{cc}1& 0^T\\ 0_{N-1}& {\tilde{D}}_{\mathrm{\Φ}}\end{array}\right]}^T---\left(9\right)$

$=\mathrm{\Φ}\left(\begin{array}{c}\left[\begin{array}{cc}\tilde{X}& x_1\end{array}\right]\end{array}\right)\left[\begin{array}{cc}\tilde{A}& -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B\\ 0^T& \frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}\end{array}\right]\×U_{\mathrm{\Φ}}^{\′}{\mathrm{\Σ}}_{\mathrm{\Φ}}^{\′}{D_{\mathrm{\Φ}}^{\′}}^T{\left[\begin{array}{cc}1& 0^T\\ 0_{N-1}& {\tilde{D}}_{\mathrm{\Φ}}\end{array}\right]}^T$

则有C^F的主要因素P可以表示为

$P=\mathrm{\Φ}\left(\begin{array}{c}\left[\begin{array}{cc}{x}_1& \tilde{X}\end{array}\right]\end{array}\right)\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B& \tilde{A}\end{array}\right]U_{\mathrm{\Φ}}^{\′}$

$=\mathrm{\Φ}\left(X\right)\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B& \tilde{A}\end{array}\right]U_{\mathrm{\Φ}}^{\′}---\left(10\right)$

又因为P＝Φ(X)A，所以得到式(11)

$A=\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B& \tilde{A}\end{array}\right]U_{\mathrm{\Φ}}^{\′}---\left(11\right)$

根据式(11)，可由下式计算得到

$\left[\begin{array}{cc}\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}& 0^T\\ -\frac{1}{h_{\mathrm{\Φ}}}\sqrt{\frac{N-1}{N(N-2)}}B& \tilde{A}\end{array}\right]=A{\left(U_{\mathrm{\Φ}}^{\′}\right)}^{-1}$

利用步骤2建立初始监测模型得到的主要因素P，计算青霉素发酵过程在特征空间中的得分向量t∈R^r为

$t=P^T[\mathrm{\Φ}\left(x_{\mathrm{new}}\right)-m_{\mathrm{\Φ}}]=A^T\mathrm{\Φ}{\left(X\right)}^T[\mathrm{\Φ}\left(x_{\mathrm{new}}\right)-\frac{1}{N}\mathrm{\Φ}\left(X\right)1_N]$

$=A^T[k(X,x_{\mathrm{new}})-\frac{1}{N}K\left(X\right)1_N]---\left(12\right)$

其中，P＝[p₁，p₂，…，p_r]，r为保留的核主元数目，k(X，x_new)表示样本矩阵X与新样本x_new的内积向量，新样本x_new的T²和SPE统计量由式(13)和式(14)计算

$T_1^2=t^T{\mathrm{\Λ}}^{-1}t---\left(13\right)$

SPE₁＝[Φ(x_new)-m_Φ]^T(I-PP^T)[Φ(x_new)-m_Φ]    (14)

其中，Λ为主元的方差阵，T²满足F分布：

$T^2=\frac{r(N^2-1)}{N(N-r)}{F}_{r,N-r}$

其中，r是保留主元的个数，T²统计量的上限为

$T_{\mathrm{\β}}^2=\frac{r(N^2-1)}{N(N-r)}{F}_{r,N-r,\mathrm{\β}}---\left(15\right)$

其中，β为置信度，而Q统计量的分布满足χ²分布，其控制上限为

Q_β＝gχ²(h)    (16)

其中，g＝ρ²/2μ，h＝2μ²/ρ²，μ与ρ²分别对应样本Q统计量的均值和方差。

3.根据权利要求1所述的基于递归核主元分析的青霉素发酵过程故障监测方法，其特征在于：步骤3所述的采用递归核主元分析方法更新步骤2的青霉素发酵过程初始的监测模型，并计算更新后青霉素发酵过程故障的主要因素

方法如下：

x_new为一个青霉素发酵过程中新样本并且可用，Φ(x_new)为新样本x_new在特征空间中的投影， $\mathrm{\Φ}\left(X_{\mathrm{new}}\right)=\mathrm{\Φ}\left(\begin{array}{c}\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]\end{array}\right)$ 为更新后青霉素发酵过程的样本矩阵，则Φ(X_new)的均值向量

和协方差阵

分别为

${\hat{m}}_{\mathrm{\Φ}}=\frac{1}{N}\mathrm{\Φ}\left(\begin{array}{c}\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]\end{array}\right)1_N=\frac{N-1}{N}{\tilde{m}}_{\mathrm{\Φ}}+\frac{1}{N}\mathrm{\Φ}\left(x_{\mathrm{new}}\right)---\left(17\right)$

${\hat{C}}^F=\frac{1}{N-1}\stackrel{\‾}{\mathrm{\Φ}}\left(\begin{array}{c}\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]\end{array}\right)\stackrel{\‾}{\mathrm{\Φ}}{\left(\begin{array}{c}\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]\end{array}\right)}^T$

$=\frac{N-2}{N-1}\left[\begin{array}{cc}\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(x_{\mathrm{new}}\right)-{\tilde{m}}_{\mathrm{\Φ}})& \sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)\end{array}\right]\×---\left(18\right)$

${\left[\begin{array}{cc}\sqrt{\frac{N-1}{N(N-2)}}(\mathrm{\Φ}\left(n_{\mathrm{new}}\right)-{\tilde{m}}_{\mathrm{\Φ}})& \sqrt{\frac{1}{N-2}}\stackrel{\‾}{\mathrm{\Φ}}\left(\tilde{X}\right)\end{array}\right]}^T$

由式(2)-(9)得到：

$V_{\mathrm{\Φ}}^{\′}=\left[\begin{array}{cc}{\widetilde{\mathrm{\Σ}}}_{\mathrm{\Φ}}& \widetilde{\mathrm{\Λ}}{\tilde{A}}^T\sqrt{\frac{N-1}{N(N-2)}}(k(\tilde{X},x_{\mathrm{new}})-\frac{1}{N-1}K\left(\tilde{X}\right)1_{N-1})\\ 0^T& h_{\mathrm{\Φ}}^{\′}\end{array}\right]$

对V′_Φ作奇异值分解得到 $V_{\mathrm{\Φ}}^{\′}=U_{\mathrm{\Φ}}^{\′\′}{\mathrm{\Σ}}_{\mathrm{\Φ}}^{\′\′}{D_{\mathrm{\Φ}}^{\′\′}}^T$

进而可以得到

的主要因素

和特征值阵

$\hat{P}=\mathrm{\Φ}\left(\begin{array}{c}\left[\begin{array}{cc}\tilde{X}& x_{\mathrm{new}}\end{array}\right]\end{array}\right)\left[\begin{array}{cc}\tilde{A}& -\frac{1}{h_{\mathrm{\Φ}}^{\′}}\sqrt{\frac{N-1}{N(N-2)}}{B}^{\′}\\ 0^T& \frac{1}{h_{\mathrm{\Φ}}^{\′}}\sqrt{\frac{N-1}{N(N-2)}}\end{array}\right]U_{\mathrm{\Φ}}^{\′\′}=\mathrm{\Φ}\left(X_{\mathrm{new}}\right)\hat{A}---\left(19\right)$

$\widehat{\mathrm{\Λ}}=\frac{N-2}{N-1}{{\mathrm{\Σ}}_{\mathrm{\Φ}}^{\′\′}}^2---\left(20\right)$

其中：青霉素发酵过程更新监测模型主要因素的调节因子：

$h_{\mathrm{\Φ}}^{\′}=\frac{N-1}{N(N-2)}\sqrt{1-2B^{\′T}k\left(\widetilde{X,}{x}_{\mathrm{new}}\right)+B^{\′T}K\left(\tilde{X}\right)B^{\′}}---\left(21\right)$

青霉素发酵过程更新监测模型主要因素的修正矩阵：

$B^{\′}=\frac{1}{N-1}{1}_{N-1}+\tilde{A}\widetilde{\mathrm{\Λ}}{\tilde{A}}^T(k\left(\widetilde{X,}{x}_{\mathrm{new}}\right)-\frac{1}{N-1}K\left(\tilde{X}\right)1_{N-1})---\left(22\right)$

表示青霉素发酵过程中间矩阵和新采样x_new的内积。

Images