CN110866643B - 基于最大二次互信息准则回归的发酵过程质量变量预测方法 - Google Patents

Abstract

本发明公开了一种基于最大二次互信息准则回归的发酵过程质量变量预测方法。实际生产中，反映最终产品质量的关键质量变量通常难以在线测量。目前常采用的离线测量法会导致滞后和精度不足等问题，影响产量及质量的一致性。本发明针对发酵生产数据强非线性、非高斯性的特点，提出一种基于最大二次互信息准则的回归方法，实现对发酵过程关键质量变量的预测。与MPLS等基于二阶统计量的回归方法相比，本发明使用高阶统计量进行过程变量与关键质量变量之间的回归，挖掘过程变量与关键质量变量之间的非线性依赖关系，且本发明无需假设数据服从高斯分布。相关实验表明，与MPLS方法相比，本发明预测效果更优。

Description

基于最大二次互信息准则回归的发酵过程质量变量预测方法

技术领域

本发明涉及基于数据驱动的回归预测方法领域，特别是涉及一种基于针对发酵生产数据的最大二次互信息准则回归方法。

背景技术

间歇过程是现代工业生产的重要生产方式。间歇过程广泛存在于生物制药、食品加工、化工、半导体生产等领域，最大的特征是小批量、高附加值、多规格和高品质。发酵过程是典型的间歇生产过程。我国生物发酵行业产业规模不断扩大，已成为我国战略性新兴产业中的重要组成部分。近年来增强生物发酵产业的自主创新能力，推动高新技术改造传统制造技术，已经成为生物发酵产业的研究热点。

在实际发酵生产过程中，一些可以反映最终产品质量的关键质量变量通常难以在线测量，如菌体浓度、产物浓度、葡萄糖浓度等。目前，这些变量大都只能离线测量，即在生产现场取样后拿到实验室进行分析测量。这样会导致滞后和精度不足等问题，影响发酵过程的产量及质量的一致性。随着传感器技术、智能仪表及计算机技术的迅速发展，生产中积累了大量的过程数据，如温度、压力及通风速率等。通过使用多元统计分析等数据驱动的方法，挖掘可测量的过程数据所包含的有用信息，实现对难以测量的变量的预测，已成为经济可靠、响应迅速的解决方式。目前，以偏最小二乘回归(Multiway Partial least squares,MPLS)及主成分回归(Principal component regression,PCR)为代表的多元统计回归模型在间歇过程的关键质量变量预测中得到广泛应用。

发酵过程的动力学模型呈现高度的非线性，因此难以建立精确的数学模型。发酵过程呈现出较强的时变性，过程的动力学特征随着发酵时间或批次不同而不断变化。此外，发酵生产数据多为高斯和非高斯的混合分布。然而，MPLS及PCR都是基于二阶统计量的方法，通过最大化方差、最大化相关性来提取特征，同时需要假设数据服从高斯分布。互信息(Mutual information)可以衡量数据分布之前的非线性相互依赖性，通过数据的概率分布提取高阶统计量。互信息在特征选择和特征提取领域得到广泛应用。为了便于计算，一些学者提出二次互信息(Quadratic mutual information,QMI)的概念。QMI将二次Renyi熵与Parzen窗密度估计方法巧妙结合，提供了有效的、便捷的计算数据集之间非线性依赖关系的方法。

发明内容

本发明针对发酵过程数据强非线性、非高斯性的特点，提出一种基于最大二次互信息准则的回归方法，实现对发酵过程关键质量变量的预测。本发明无需假设数据服从高斯分布，通过二次互信息提取过程变量与关键质量变量之间的非线性关系，提高预测精度。

本发明的总体设计方案为：

首先，将高维的原始过程数据(输入空间)进行线性映射到低维的特征空间。如图1所示，通过定义目标函数，使得特征与关键质量变量(输出)之间的二次互信息最大，同时使输入空间的Renyi熵也尽量保持最大。然后，使用梯度下降法求取最优线性变换。最后，建立原始过程数据与关键质量变量之间的线性回归模型，对新样本时刻的关键质量变量进行预测。

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

1、基于最大二次互信息准则回归的发酵过程质量变量预测，其特征在于：包括以下步骤：

A.离线建模阶段：

1)三维数据展开；正常发酵生产中采集到的三维过程数据可表示为X_o(I×M×K)，其中I为批次数，M为过程变量个数，K为采样时刻数，首先对其进行变量展开，经过变量展开后得到二维矩阵X_v(N×M),其中N＝KI；对于采集到的三维关键质量变量矩阵进行相同的处理，得到Y_v(N×L)，其中L为关键质量变量的个数；最后对X_v和Y_v进行转置运算，得到X＝X_v ^T＝[x₁,…,x_N]及Y＝Y_v ^T＝[y₁,…,y_N]；

2)数据标准化；对X和Y进行标准化处理，求取每一行的均值和标准差；X的均值为：

其中

X的标准差为：

X_std＝[x_s1,...,x_si,x_sM]^T (2)

其中

对Y进行标准化处理，得到Y的均值为：

其中

Y的标准差为：

Y_std＝[y_s1,...,y_si,y_sL]^T(4)

其中

3)依据Silverman’s rule计算使用Parzen窗进行密度估计时所需的核函数，包括σ_x,σ_t与σ_y；首先计算标准化后X的协方差矩阵Σ，令

其中σ_ii为协方差矩阵Σ对角线上元素；依据Silverman’srule：

/>

其中r为对X进行线性变化后特征空间的维数；

4)目标函数J(W)的定义及计算；对X进行线性变换T＝W^TX，其中W是M×r(r<M)维变换矩阵，W的每一列为一个投影轴；T为变换后的特征空间,t_i(i＝1,…,N)为x_i在特征空间的特征向量；即

目标函数定义为：

J(W)＝-H_R2(T)-α·QMI_ED(T,Y) (9)

其中H_R2(T)是T的二次Renyi熵，QMI_ED(T,Y)是X与Y之间的基于欧式距离的二次互信息，α是权重系数，J(W)取得最小值,即W＝W*时，T与Y之间的二次互信息最大，同时使T的二次Renyi熵也尽量保持最大；

使用Parzen窗对H_R2(T)和QMI_ED(T,Y)进行计算；首先令高斯核函数记为：

然后使用Parzen窗对H_R2(T)和QMI_ED(T,Y)分别进行估计，计算方法如下：

QMI_ED(T,Y)＝V_E＝V_J+V_M-2V_C

/>

5)使用梯度下降法求取W*；随机初始化W，使用3)4)中相关公式计算此时J(W)；使用

对W值进行更新，再计算J(W)，其中，

以此反复进行迭代计算，直到达到终止条件为止，其中λ为学习速率；终止条件可以设置为达到最大迭代次数或J(W)收敛至某精度；

6)建立回归模型；由5)求得W*,此时T＝W*^TX；则T与Y之间的回归系数为：θ＝(TT^T)^-1TY^T；

B.在线预测阶段：

1)样本预处理；对于采集到的新样本x_new，对其进行正则化处理：

其中x_new,i表示x_new的第i个分量；同理X_aver,i表示X_aver的第i个分量,X_std,i表示X_std的第i个分量；

2)使用回归模型进行预测

y_new＝θ^TW*^Tx_new (14)

3)将y_new映射回原始输出空间，求取最终预测结果

y_p,i＝y_new,i·Y_std,i+Y_aver,i (i＝1,…,L) (15)

y_p,i表示y_p的第i个分量，Y_aver,i表示Y_aver的第i个分量,Y_std,i表示Y_std的第i个分量。

有益效果

本发明针对发酵过程数据强非线性、非高斯性的特点，提出一种基于最大二次互信息准则的回归方法，实现对发酵过程中难以控制的关键质量变量的预测，提高预测精度。与现有的技术PLS或PCR相比，本发明无需假设数据服从高斯分布。而且，本发明依据高阶统计量，即二次互信息，提取过程变量与关键质量变量之间的非线性依赖关系，与仅考虑二阶统计量、提取线性相关性的MPLS或PCR方法相比，具有明显优势。

附图说明

图1求取最优线性变换示意图；

图2三维数据变量展开示意图；

图3本发明与MPLS方法对某批次OD值的预测结果对比图；

具体实施方式

现代生物制药领域中，通常将目的基因导入菌体内，形成基因工程菌，经过培养与发酵使其表达，以生成所需的药用蛋白。大肠杆菌是常用的基因工程菌之一。实际药物制备过程中，大肠杆菌常被基因改良并进行发酵生产，以制备重组人白介素-2(IL-2)。IL-2是一种重要的药用蛋白，广泛应用于恶性肿瘤的治疗。大肠杆菌的发酵过程是典型的间歇生产过程，包含一系列复杂的生物化学反应，数据具有高度的非线性及高斯性。实际生产数据验证更突显本发明的研究意义和效果，因此本发明以实际制备IL-2的生产数据为验证对象。实际生产中，常通过离线测量光密度(optical density,OD)间接反映菌体浓度，判断菌体生长状态，OD是关键质量变量。

实验中，共计选择8个可测量的生产过程变量(浓度罐压、温度、酸碱度、搅拌速率、补碳、补氮、通风速率、溶解氧、)与关键质量变量OD。发酵时长约为6.5小时，采样间隔为10分钟，共计39个采样时刻。本发明共收集到28批次正常生产数据。随机选取其中的20批次进行离线回归建模，剩余8批次用于测试实验。

将本发明方法应用的上述大肠杆菌发酵制备IL-2的实际生产数据集中。开发语言为Matlab，开发平台为Matlab R2014a。整个实施过程包括离线建模与在线预测两部分。

A.离线建模阶段：

1)生产数据预处理。正常发酵生产中采集到的三维数据可表示为X_o(I×M×K)，其中批次总数I＝20，过程变量个数M＝8，采样点总数K＝39，首先对其进行变量展开，如图2所示，得到二维矩阵X_v(780×8),其中样本总数N＝KI＝780。对质量变量OD矩阵进行相同的处理，得到Y_v(780×1)，其中L＝1。最后X_v和Y_v进行转置运算，得到X_8×780＝X_v ^T＝[x₁,…,x_i,…,x₇₈₀]及Y_1×780＝Y_v ^T＝[y₁,…,y_i,…,y₇₈₀]。

2)数据标准化及相关统计量计算。对X和Y进行标准化处理，求取每一行的均值和标准差。使用Matlab中“mean”函数求取X的均值X_aver(1×8)，使用“std”函数求取X的标准差X_std(1×8)，同理，对Y进行标准化处理，得到Y的均值Y_aver及标准差Y_std。计算σ_x,σ_t与σ_y。首先使用Matlab中“cov”函数计算标准化后X的协方差矩阵Σ，求取Σ对角线上元素的均值

依据Silverman’s rule分别计算得到：

令r＝1,计算得到

3)使用梯度下降法求取最优线性变换。对X进行线性变换T＝W^TX，r＝1时，对应的W是8×1维变换矩阵。随机初始化W。计算此时特征向量t_i＝W^Tx_i,即

T＝W^TX

[t₁ t_i … t₇₈₀]＝W^T[x₁ x_i… x₇₈₀]

计算目标函数J(W)。首先使用Parzen窗对H_R2(T)和QMI_ED(T,Y)分别进行估计，计算方法如下：

/>

QMI_ED(T,Y)＝V_E＝V_J+V_M-2V_C

计算目标函数J(W),权重系数α＝80。使用

对W值进行更新，再计算J(W)。以此反复进行迭代计算，直到达到终止条件为止。其中学习速率λ＝0.5。终止条件设置为J(W)收敛至精度0.001。

4)建立回归模型。由5)求得W*,此时T＝W*^TX。则T与Y之间的回归系数为：θ＝(TT^T)^-1TY^T。

B.在线预测阶段：

1)新样本预处理。8批次的测试数据集中的每个批次数据分别进行单独处理。针对某一个测试批次的数据，可表示为X_test(M×K)。其中，过程变量个数M＝8，采样点总数K＝39。X_test的每一列表示测试批次在某采样点所有变量的取值，即X_test的每一列为一个新样本。对于采集到的每个新样本x_new(8×1)，对其进行正则化处理,使用Matlab中的点除函数，计算如下：

x_new＝(x_new-X_aver)./X_std

2)使用回归模型进行预测。y_new＝θ^TW*^Tx_new

3)将y_new映射回原始输出空间，求取最终预测结果。使用Matlab中的点乘函数。得到最终预测结果:y_p＝y_new.*Y_std+Y_aver

使用以上方法对测试批次的所有采样时刻的样本进行预测。进一步对总计8批次的测试数据逐个进行预测，得到最终预测结果。

表1为本发明与传统MPLS方法对8个测试批次预测结果的预测均方根误差指标RMSE(Root mean square error,RMSE)值。从表中数据看出，本发明的预测精度要优于传统的MPLS方法。图3为某个测试批次下两种方法对OD值的预测结果。可以看出，在采样时刻21附近及批次结束时，本发明预测误差明显小于MPLS,具有更好的预测效果。

表1

测试批次	1	2	3	4	5	6	7	8
									MPLS	10.5515	10.3351	10.305	11.8488	11.0787	10.0776	9.3867	8.3963
本发明	8.9149	7.8731	9.5526	6.6075	9.9596	8.8506	6.7318	7.3542

Claims (1)

1.基于最大二次互信息准则回归的发酵过程质量变量预测方法，其特征在于：包括以下步骤：

A.离线建模阶段：

1)三维数据展开；正常发酵生产中采集到的三维过程数据可表示为X_o(I×M×K)，其中I为批次数，M为过程变量个数，K为采样时刻数，首先对其进行变量展开，经过变量展开后得到二维矩阵X_v(N×M),其中N＝KI；对于采集到的三维关键质量变量矩阵进行相同的处理，得到Y_v(N×L)，其中L为关键质量变量的个数；最后对X_v和Y_v进行转置运算，得到X＝X_v ^T＝[x₁,…,x_N]及Y＝Y_v ^T＝[y₁,…,y_N]；

2)数据标准化；对X和Y进行标准化处理，求取每一行的均值和标准差；X的均值为：

其中

X的标准差为：

X_std＝[x_s1,...,x_si,x_sM]^T (2)

其中

对Y进行标准化处理，得到Y的均值为：

其中

Y的标准差为：

Y_std＝[y_s1,...,y_si,y_sL]^T (4)

其中

3)依据Silverman’s rule计算使用Parzen窗进行密度估计时所需的核函数，包括σ_x,σ_t与σ_y；首先计算标准化后X的协方差矩阵Σ，令

其中σ_ii为协方差矩阵Σ对角线上元素；依据Silverman’s rule：

其中r为对X进行线性变化后特征空间的维数；

4)目标函数J(W)的定义及计算；对X进行线性变换T＝W^TX，其中W是M×r r<M维变换矩阵，W的每一列为一个投影轴；T为变换后的特征空间,t_i i＝1,…,N为x_i在特征空间的特征向量；即

目标函数定义为：

J(W)＝-H_R2(T)-α·QMI_ED(T,Y) (9)

其中H_R2(T)是T的二次Renyi熵，QMI_ED(T,Y)是X与Y之间的基于欧式距离的二次互信息，α是权重系数，J(W)取得最小值,即W＝W*时，T与Y之间的二次互信息最大，同时使T的二次Renyi熵也尽量保持最大；

使用Parzen窗对H_R2(T)和QMI_ED(T,Y)进行计算；首先令高斯核函数记为：

然后使用Parzen窗对H_R2(T)和QMI_ED(T,Y)分别进行估计，计算方法如下：

QMI_ED(T,Y)＝V_E＝V_J+V_M-2V_C

5)使用梯度下降法求取W*；随机初始化W，使用3)4)中相关公式计算此时J(W)；使用

对W值进行更新，再计算J(W)，其中，

以此反复进行迭代计算，直到达到终止条件为止，其中λ为学习速率；终止条件设置为达到最大迭代次数或J(W)收敛至某精度；

6)建立回归模型；由5)求得W*,此时T＝W*^TX；则T与Y之间的回归系数为：θ＝(TT^T)^-1TY^T；

B.在线预测阶段：

1)样本预处理；对于采集到的新样本x_new，对其进行正则化处理：

其中x_new,i表示x_new的第i个分量；同理X_aver,i表示X_aver的第i个分量,X_std,i表示X_std的第i个分量；

2)使用回归模型进行预测

y_new＝θ^TW*^Tx_new (14)

3)将y_new映射回原始输出空间，求取最终预测结果

y_p,i＝y_new,i·Y_std,i+Y_aver,i i＝1,…,L(15)

y_p,i表示y_p的第i个分量，Y_aver,i表示Y_aver的第i个分量,Y_std,i表示Y_std的第i个分量。

Images