CN113031553A - 一种基于自动聚类结合偏最小二乘的间歇过程质量预测方法 - Google Patents

Abstract

本发明公开了一种基于自动聚类结合偏最小二乘的间歇过程质量预测方法。针对间歇过程的多阶段性特征，目前已有的阶段划分方法很少考虑质量相关变量对阶段划分结果的影响。本发明在划分阶段前利用典型相关分析(canonical correlation analysis，CCA)对间歇过程数据进行特征选择，在保证其过程变量以及质量相关变量之间相关关系最大时找到其最优的线性表示。该过程不仅可以实现数据降维，同时考虑质量相关变量对划分结果的影响。最终，在DBSCAN划分阶段内建立基于MPLS的质量预测模型。将该算法在青霉素发酵仿真实验平台进行了实验验证，实验结果证明了本方法的可行性和有效性。

Description

一种基于自动聚类结合偏最小二乘的间歇过程质量预测方法

技术领域

本发明涉及基于数据驱动的多元统计过程监控(multivariate statisticalprocess monitoring,MSPM)技术领域，特别是针对间歇过程数据的多阶段特性，提出一种基于自动聚类结合偏最小二乘的间歇过程质量预测方法。

背景技术

间歇过程又称分批过程，是精细化工、生物制药以及食品饮料等行业中最为主要的生产方式，诸如青霉素发酵，大肠杆菌制备药用蛋白质等都属于间歇过程范畴。对间歇过程进行质量预测，依据预测值进一步指导操作，控制生产过程，确保产品质量符合生产标准十分具有必要性。

间歇过程复杂多变，依靠其机理知识进行建模往往难以实现，传感器技术的飞速发展使得诸如偏最小二乘(partial least-square，PLS)、主成分分析(principalcomponent analysis，PCA)等基于数据驱动的多元统计方法被广泛应用于间歇过程的质量预测研究。其中，多向偏最小二乘算法(multiway partial least-square，MPLS)作为PLS的一种拓展，被视为间歇过程质量预测研究领域的重大突破。多阶段特性是间歇生产过程固有特性之一，不同阶段内其数据特征存在明显差异。针对间歇过程的多阶段特性，国内外学者开展了诸多研究工作，Zhao等提出了基于K-均值的间歇过程时段划分方法，该方法引入模糊隶属度作为两个相邻稳定阶段的权重系数，用相邻稳定阶段来近似表示二者之间的过渡阶段，一定程度提高了模型的精度。Gao等提出了一种基于扩展核熵负载矩阵的阶段划分方法，该方法首先利用特征提取后的主元个数进行初步划分，再将其负载矩阵进行扩展并利用模糊c均值进行二次划分。间歇过程种类繁多，诸如青霉素发酵过程、大肠杆菌制备药用蛋白质等皆属于间歇过程范畴，K-均值、模糊c均值等聚类方法在进行阶段划分之前需要人为指定簇类个数，然而对于一个陌生、复杂的间歇过程，研究人员难以得知其确切阶段个数。Ge等提出了一种基于重复性因子的多阶段在线质量预测方法；邓小刚等利用时序核熵主元关联矩阵的相似性实现阶段划分。然而，上述方法在进行阶段划分时均只关注过程变量，忽视质量相关变量对划分结果的影响。

发明内容

本发明针对间歇过程的多阶段特性，提出一种基于自动聚类结合偏最小二乘的间歇过程质量预测方法。基于密度的可剔除噪声点的聚类(density-based spatialclustering of applications with noise，DBSCAN)是一种自动聚类的方法，不仅无需事先为指定簇的个数，在机理知识未知的间歇过程领域亦可发挥作用，而且可以有效识别出噪声点，尤为适用于工况复杂的间歇过程。本发明在划分阶段前利用典型相关分析(canonical correlation analysis，CCA)对间歇过程数据进行特征选择，在保证其过程变量以及质量相关变量之间相关关系最大时找到其最优的线性表示。该过程不仅可以实现数据降维，同时考虑质量相关变量对划分结果的影响。最终，在DBSCAN划分阶段内建立基于MPLS的质量预测模型。

间歇过程数据存在固有的三维特性，假设模型建立所采用的各批次操作时间皆等长，则可得过程数据X(I×J_X×K)以及质量数据Y(I×J_Y×K),其中I为批次个数，K为每批次的采样点个数，J_X为过程数据的变量个数，J_Y则是由离线测量得到的质量相关数据的变量个数。则本发明具体步骤如下：

Step1：将X(I×J_X×K)沿批次方向展开，得到K个时间片子矩阵X_k(I×J_X)并对其进行标准化处理，即：

其中，

即为标准化后的数据，

为第k采样时刻的第j过程变量在所有批次过程中的均值，S_k,j即为其标准差。质量相关数据Y(I×J_Y×K)亦做标准化处理；

Step2：对各时间片矩阵进行基于CCA的特征提取过程，以期得到在过程数据和质量相关数据相关关系最大时的最优线性表示X′和Y′，使得二者之间的相关系数ρ(X′,Y′)最大,其中a与b即为线性系数向量，上标T则为转置含义，即:

X′＝a^TX,Y′＝b^TY (2)

由相关系数定义易知优化目标为求取在相关系数

最大时，a与b的值，即：

其中cov(X′,Y′)为X′和Y′的协方差，D(X′)以及D(Y′)则分别为二者的方差。由于X和Y均值为0且方差为1，并令协方差cov(X,Y)＝S_XY，则优化目标为：

其中，S_XX表示为矩阵X的方差，S_YY表示为矩阵Y的方差，采用与SVM类似的优化方法，固定分母，优化分子，满足a^TS_XXa＝1,b^TS_YYb＝1，则优化目标可转化为：

s.t.a^TS_XXa＝1,b^TS_YYb＝1

令线性系数向量a＝S_XX ^-1/2u,b＝S_YY ^-1/2v，并满足u^Tu＝1,v^Tv＝1则优化目标最终为：

s.t.u^Tu＝1,v^Tv＝1

为求解该式，对M＝S_XX ^-1/2S_XYS_YY矩阵做奇异值分解有：

u^TS_XX ^-1/2S_XYS_YY ^-1/2v＝u^TUΣV^Tv (6)

其中U和V分别为M矩阵的左奇异向量u和右奇异向量v组成的矩阵，Σ为其奇异值组成的对角矩阵。因此，Σ中最大奇异值即为优化目标的最大值，利用其对应的左右奇异向量u和v即可求出线性系数a＝S_XX ^-1/2u,b＝S_YY ^-1/2v，进一步由可由X′＝a^TX,Y′＝b^TY得到最优线性表示X′和Y′。将各批次间歇过程过程变量以及质量相关变量的最优线性表示组合为扩展矩阵[X′Y′]，其中第k时刻表示为P_k(I×2),k＝1,2,...,K，由于与Step3中所提到的k-dist定义冲突，因此将其记为p(i)，K个采样点所构成的数据集记为P＝{p(i)}；

Step3:设置DBSCAN算法中邻域半径阈值Eps，即为对于给定点p为中心的圆形邻域范围。

首先明确k-dist的定义是指对于给定数据集P＝{p(i)}其中i＝0,1,...,n，对于任意点p(i)，计算其到集合S＝{p(1),p(2),...,p(i-1),p(i+1),...,p(n)}中所有点之间的距离，将从小到大进行排序，得到集合D＝{d(1),d(2),...,d(k-1),d(k),d(k+1),...,d(n)}，k值可设置为4，其中d(k)为点p(i)到除自己外其余所有点之间距离第k近的距离，即k-dist，对待聚类集合中的每个点都计算其k-dist，最终得到所有的点的k-dist集合E＝{e(1),e(2),...,e(n)}，将其进行升序排列得到E′，拟合一条排序后的变化曲线图，选取其中变化最为剧烈的位置所对应的k-dist为Eps的值。

Step4：优化DBSCAN算法中最小包含点数MinPts，即Step3中的k值。本发明使用轮廓系数(Silhouette Coefficient)来进行求取，样本轮廓系数定义如下：

对于点b(i)，a(i)即为该点到与其同一类别其余样本之间的平均距离，b(i)则为该点与不同类别中样本之间的平均距离的最小值。通过调节最小包含点数MinPts参数大小，使得被分类样本所有样本点s(i)→1，并相应调节Step3中设置邻域半径阈值Eps时所选取的k值。

Step5：选择欧氏距离作为DBSCAN算法的距离度量，即对于点p(a)和p(b)所对应的A,B矩阵，其组成元素分别表示为a_ij以及b_ij，a_ij与b_ij分别对应A,B矩阵中第i行，第j列的值，其中i＝1,2,...I,j＝1,2，则点p(a)和点p(b)欧式距离定义如下：

Step6：依据Step4、Step5设置DBSCAN算法的超参数Eps以及MinPts，以进行特征提取后的样本作为输入数据，基于DBSCAN进行聚类得到最终的分类结果C＝{C₁,C₂,...,C_N}，C_n,n＝1,2,...,N即代表DBSCAN聚类后所形成的第n簇，即间歇过程的第n阶段，基于DBSCAN的聚类具体步骤如下：

6.1初始化核心对象集合Ω＝φ，最终聚类结果数N＝0，未被访问样本集合Γ＝D，簇划分C＝φ；

6.2对于K个采样点所构成的数据集记为P＝{p(i)}，i＝1,2,...,K，依据以下步骤找出其所有核心对象：

(1)依据距离度量方式，找到样本x_k的领域半径阈值内所有子样本集N_Eps(x_k)；

(2)若子样本集个数满足|N_Eps(x_k)|≥MinPts，将样本x_k加入核心对象样本集合Ω＝Ω∪{x_j}；

6.3若核心对象集合Ω＝φ，则算法结束，否则继续执行6.4；

6.4在核心对象Ω集合中，随机选取任意核心对象o，初始化当前簇的核心对象队列为Ω_cur＝{o}，初始化类别序号为N＝N+1,初始化当前簇的样本集合为C_N＝{o}，更新未被访问的样本集合Γ＝Γ-{o}；

6.5若当前簇的核心对象队列Ω_cur＝φ，则C_c聚类完毕，更新划分结果C＝{C₁,C₂,...,C_N}，更新核心对象集合Ω＝Ω-C_N，转入6.3，否则更新核心对象集合Ω＝Ω-C_N；

6.6在当前簇的核心对象队列Ω_cur中取出一个核心对象o′，找出所有在领域半径阈值Eps范围内的子样本集N_Eps(o′)，令Δ＝N_Eps(o′)∩Γ，更新当前簇的样本集合C_N＝C_N∪Δ，更新未被访问的样本集合Γ＝Γ-Δ，更新Ω_cur＝Ω_cur∪(Δ∩Ω)-o′，转入6.5；

Step7：依据划分结果在各个子阶段分别建立基于MPLS的回归预测模型，对于在线采集到的数据，依据其采样时间判断其所属阶段，利用该阶段所建立模型进行质量预测；

对于第n阶段的矩阵，MPLS将标准化后的数据矩阵X_n和Y_n进行如下分解，其中X_n和Y_n表示截取X和Y中第n阶段的子矩阵：

X_n＝TP^T+E (9)

Y_n＝UQ^T+F (10)

其中，T和U分别为X_c和Y_c的主元得分矩阵，P和Q为负载矩阵，E和F为负载矩阵，则第n阶段的最终回归预测模型可表示为：

Y_n＝X_nβ+E_Y (11)

其中β即为质量预测的系数矩阵，E_Y则为残差矩阵。

有益效果

本发明在对间歇过程进行阶段划分时，不仅无需人为指定簇类个数，在复杂陌生的间歇过程领域亦适用，而且同时考虑了过程变量以及质量相关变量对划分结果的影响，有效提高阶段划分的结果精度，进而进一步提高质量预测结果精度，对间歇过程过程的质量预测存在重要意义。

附图说明

图1所示为本发明流程图示；

图2所示为阶段划分结果图示；

图3所示为16个测试批次菌体浓度的预测误差；

图4所示为16个测试批次产物浓度的预测误差；

图5所示为某一测试批次产物浓度预测效果图示；

图6所示为某一测试批次菌体浓度预测效果图示；

具体实施方式

青霉素发酵过程是工业上一种典型间歇过程，其生产过程包括：快速消耗底物促进菌体生长时期以及合成青霉素时期，两者分别耗时约为45h以及355h。Pensim仿真平台由伊利诺伊理工学院(Illinois Institute of Technology，IIT)研发，是国内网上诸多学者用于仿真青霉素发酵过程的实验平台，其可以通过改变青霉素生产过程的微生物浓度、CO2浓度、氧浓度、pH值等变量来模拟不同操作条件。设定采样时间间隔为1h，如表1所示，选取10个过程变量以及2个质量相关变量，利用Pensim2.0仿真平台进行56批的工况仿真，得到三维数据X(56×10×400)以及Y(56×2×400)，其中40批作为训练数据，16批作为测试数据。

表1青霉素发酵过程变量

Table 1 Process variables of Penicillin fermentation

基于以上描述，按照发明内容，将具体过程在pycharm中实现如下：

Step1：依据发明内容Step1将训练数据X(40×10×400)以及Y(40×2×400)沿批次方向展开，得到时间片矩阵X_i(10×400)以及Y_i(10×400)，其中i＝1,2,...,I，并对其进行标准化处理，标准化的数据记为X以及Y；

Step2：对各时间片矩阵进行基于CCA的特征提取过程，以期得到在过程数据和质量相关数据相关关系最大时的最优线性表示X′＝a^TX,Y′＝b^TY，使得二者之间的相关系数ρ(X′,Y′)最大，将各批次间歇过程过程变量以及质量相关变量的最优线性表示组合为扩展矩阵[X′Y′]，其中第k时刻表示为P_k(400×2),k＝1,2,...,K，由于与Step3中所提到的定义冲突，因此将其记为p(i)，K个采样点所构成的数据集记为P＝{p(i)}；

Step3：依据发明内容Step3所述，设置DBSCAN算法中参数邻域半径阈值Eps；

Step4：依据发明内容Step3所述，设置DBSCAN算法中参数最小包含点数MinPts；

Step5：依据发明内容Step3所述，选择欧氏距离作为DBSCAN算法的距离度量；

Step6：依据Step4、Step5设置DBSCAN算法的超参数Eps以及MinPts，以进行特征提取后的样本作为输入P＝{p(i)}数据，基于DBSCAN进行聚类得到阶段划分结果C＝{C₁,C₂,...,C_k}。

Step7：依据划分结果在各个子阶段分别建立基于MPLS的回归预测模型，对于在线采集到的数据，依据其采样时间判断其所属阶段，利用该阶段所建立模型进行质量预测；

Step8:为验证本发明所提出质量预测方法的有效性，本文与基于传统MPLS方法进行质量预测进行对比。

判断模型预测精度所采用的标准为均方根误差(Root Mean Square Error,RMSE)，即

其中，K为当前预测批次的采样点个数，

表示第k采样点的预测值，y_k则表示第k采样点的实际值。需要注意的是，间歇过程新批次的整体预测效果计算需要等到该批次采集全部完成后进行。

由图示可看出，青霉素生产过程明显的被分成不同阶段。在与不考虑质量变量信息的阶段划分方法的对比中可以看出，考虑质量变量后的阶段划分方法在质量预测的精度和稳定性方面都有所提高。

Claims (1)

1.一种基于自动聚类结合偏最小二乘的间歇过程质量预测方法，其特征在于：

间歇过程数据存在固有的三维特性，模型建立所采用的各批次操作时间皆等长，得到过程数据X(I×J_X×K)以及质量数据Y(I×J_Y×K)，其中I为批次个数，K为每批次的采样点个数，J_X为过程数据的变量个数，J_Y则是由离线测量得到的质量相关数据的变量个数；

具体步骤如下：

Step1：将X(I×J_X×K)沿批次方向展开，得到K个时间片子矩阵X_k(I×J_X)并对其进行标准化处理，即：

其中，

即为标准化后的数据，

为第k采样时刻的第j过程变量在所有批次过程中的均值，S_k,j即为其标准差。质量相关数据Y(I×J_Y×K)亦做标准化处理；

Step2：对各时间片矩阵进行基于CCA的特征提取过程，以期得到在过程数据和质量相关数据相关关系最大时的最优线性表示X′和Y′，将各批次间歇过程过程变量以及质量相关变量的最优线性表示组合为扩展矩阵[X′Y′]，其中第k时刻表示为P_k(I×2),k＝1,2,...,K，由于与Step3中所提到的k-dist定义冲突，因此将其记为p(i)，K个采样点所构成的数据集记为P＝{p(i)}；

Step3:设置DBSCAN算法中邻域半径阈值Eps，即为对于给定点p为中心的圆形邻域范围，首先明确k-dist的定义是指对于给定数据集P＝{p(i)}其中i＝0,1,...,n，对于任意点p(i)，计算其到集合S＝{p(1),p(2),…,p(i-1),p(i+1),...,p(n)}中所有点之间的距离，将从小到大进行排序，得到集合＝{d(1),d(2),...,d(k-1),d(k),d(k+1),...,d(n)}，k值可设置为4，其中d(k)为点p(i)到除自己外其余所有点之间距离第k近的距离，即k-dist，对待聚类集合中的每个点都计算其k-dist，最终得到所有的点的k-dist集合E＝{e(1),e(2),...,e(n)}，将其进行升序排列得到E′，拟合一条排序后的变化曲线图，选取其中变化最为剧烈的位置所对应的k-dist为Eps的值；

Step4：优化DBSCAN算法中最小包含点数MinPts，即Step3中的k值，使用轮廓系数(Silhouette Coefficient)来进行求取，样本轮廓系数定义如下：

对于点b(i)，a(i)即为该点到与其同一类别其余样本之间的平均距离，b(i)则为该点与不同类别中样本之间的平均距离的最小值。通过调节最小包含点数MinPts参数大小，使得被分类样本所有样本点s(i)→1，并相应调节Step3中设置邻域半径阈值Eps时所选取的k值；

Step5：选择欧氏距离作为DBSCAN算法的距离度量，即对于点p(a)和p(b)所对应的A,B矩阵，其组成元素分别表示为a_ij以及b_ij，a_ij与b_ij分别对应A,B矩阵中第i行，第j列的值，其中i＝1,2,...I，j＝1,2，则点p(a)和点p(b)欧式距离定义如下：

Step6：依据Step4、Step5设置DBSCAN算法的超参数Eps以及MinPts，以进行特征提取后的样本作为输入数据，基于DBSCAN进行聚类得到最终的分类结果C＝{C₁,C₂,...,C_N}，C_n,n＝1,2,...,N即代表DBSCAN聚类后所形成的第n簇，即间歇过程的第n阶段；

Step7：依据划分结果在各个子阶段分别建立基于MPLS的回归预测模型，对于在线采集到的数据，依据其采样时间判断其所属阶段，利用该阶段所建立模型进行质量预测；

对于第n阶段的矩阵，MPLS将标准化后的数据矩阵X_n和Y_n进行如下分解，其中X_n和Y_n表示截取X和Y中第n阶段的子矩阵：

X_n＝TP^T+E (4)

Y_n＝UQ^T+F (5)

其中，T和U分别为X_c和Y_c的主元得分矩阵，P和Q为负载矩阵，E和F为负载矩阵，则第n阶段的最终回归预测模型可表示为：

Y_n＝X_nβ+E_Y (6)

其中β即为质量预测的系数矩阵，E_Y则为残差矩阵。

Images