基于核二次互信息的發(fā)酵過程質(zhì)量預(yù)測模型

李 征，王 普，高學(xué)金+，齊詠生，高慧慧

(1.北京工業(yè)大學(xué) 信息學(xué)部，北京 100124；2.北京工業(yè)大學(xué) 數(shù)字社區(qū)教育部工程研究中心，北京 100124；3.北京工業(yè)大學(xué) 城市軌道交通北京實驗室，北京 100124；4.北京工業(yè)大學(xué) 計算智能與智能系統(tǒng)北京市重點實驗室，北京 100124；5.內(nèi)蒙古工業(yè)大學(xué) 電力學(xué)院，內(nèi)蒙古 呼和浩特 010051)

0 引 言

發(fā)酵過程廣泛應(yīng)用于生物制藥領(lǐng)域。發(fā)酵過程中的質(zhì)量變量難以在線預(yù)測，通常在批次結(jié)束后取適量發(fā)酵液帶到實驗室離線分析，具有時間延遲，不利于產(chǎn)品質(zhì)量閉環(huán)控制。過程工業(yè)中常利用傳感器讀取的過程變量預(yù)測質(zhì)量變量[1-3]。發(fā)酵過程動力學(xué)模型為強(qiáng)非線性，應(yīng)用線性方法[4-6]難以保證預(yù)測精度。文獻(xiàn)[7]應(yīng)用核偏最小二乘法對關(guān)鍵變量進(jìn)行預(yù)測，通過核映射將被測變量與可測變量間的非線性關(guān)系近似線性化。文獻(xiàn)[8]使用最小二乘支持向量回歸預(yù)測青霉素發(fā)酵中的質(zhì)量變量，解決非線性問題。文獻(xiàn)[9]提出約化雙核偏最小二乘法，將核映射后的高維特征逆映射至質(zhì)量空間，實現(xiàn)發(fā)酵過程中質(zhì)量變量的預(yù)測。文獻(xiàn)[10]建立核熵成分分析所提取的過程特征與質(zhì)量數(shù)據(jù)之間的回歸模型進(jìn)行質(zhì)量預(yù)測。然而，以上方法均基于相關(guān)性提取特征，不足以描述變量之間的非線性關(guān)系。文獻(xiàn)[11-14]使用二次互信息衡量變量之間非線性關(guān)系提取特征，并解決了互信息計算量較大的難題。文獻(xiàn)[15,16]在分類問題中通過最大化特征與類別標(biāo)簽之間的二次互信息，實現(xiàn)有監(jiān)督特征提取。

為解決發(fā)酵過程質(zhì)量預(yù)測中數(shù)據(jù)強(qiáng)非線性問題，提出一種基于核二次互信息的質(zhì)量預(yù)測模型，依據(jù)信息論準(zhǔn)則有監(jiān)督地提取過程數(shù)據(jù)的特征，并建立回歸模型對質(zhì)量變量進(jìn)行在線預(yù)測。

1 相關(guān)理論

1.1 發(fā)酵過程工藝

發(fā)酵過程工藝如圖1所示。通過基因工程等方式改變微生物遺傳特性，形成生產(chǎn)用的菌種。使用搖瓶對菌種進(jìn)行擴(kuò)大培養(yǎng)，將其接種至滅菌的發(fā)酵罐中進(jìn)行生長培養(yǎng)。發(fā)酵罐中富含菌體生長所必需的營養(yǎng)原料。發(fā)酵過程包括菌種適應(yīng)期、菌體快速生長期及蛋白表達(dá)期。微生物在合適的發(fā)酵條件下會合成代謝產(chǎn)物，某些代謝產(chǎn)物經(jīng)過分離純化后得到產(chǎn)品——目的蛋白。例如，青霉菌在合適發(fā)酵條件下可代謝產(chǎn)生廣譜抗生素青霉素；大腸桿菌是現(xiàn)代生物制藥領(lǐng)域中常用的基因工程菌之一，大腸桿菌在工業(yè)生產(chǎn)中常被基因改良并進(jìn)行發(fā)酵生產(chǎn)，以制備抗腫瘤藥物重組人白介素-2(Interleukin-2, IL-2)。發(fā)酵過程包含微生物代謝生長，過程機(jī)理復(fù)雜，難以建立精確的數(shù)學(xué)模型。發(fā)酵產(chǎn)品質(zhì)量極易受到溫度、溶氧等環(huán)境條件及滅菌、接種等人為操作的影響，波動較大。對菌體數(shù)量、菌體濃度、目的蛋白等質(zhì)量變量進(jìn)行在線預(yù)測，判斷生產(chǎn)狀態(tài)，具有十分重要的意義。

圖1 發(fā)酵過程工藝

1.2 Renyi二次熵及其估計

信息論創(chuàng)始人香農(nóng)定義了著名的香農(nóng)熵，實現(xiàn)對隨機(jī)變量所包含平均信息量的量化。熵值越大，代表隨機(jī)變量的不確定性越大，所包含的信息量就越大。Renyi對熵進(jìn)行推廣，Renyi的a階熵定義如下

(1)

其中，X是隨機(jī)變量，f(x)是其概率密度函數(shù)。在a=2時，得到Renyi二次熵，表達(dá)式如下

(2)

當(dāng)a→1時，Renyi熵近似為香農(nóng)熵。Renyi二次熵在計算量上較香農(nóng)熵具有明顯優(yōu)勢。由Renyi二次熵定義式可知，Renyi二次熵中的log運算符在積分號外面，因此計算起來相對簡便。因此，基于Renyi二次熵定義的二次互信息在計算上較經(jīng)典的基于香農(nóng)熵意義的互信息具有明顯優(yōu)勢。

(3)

(4)

f(x)是產(chǎn)生數(shù)據(jù)集X的概率密度函數(shù)。 p(u,Σ) 是對f(x)進(jìn)行Parzen窗估計時選取的高斯核函數(shù)。為便于計算，令Σ=σ2I為其協(xié)方差矩陣，以此確定窗寬大小。式(3) 中利用高斯核函數(shù)的重要性質(zhì)[15]，將繁雜的積分運算轉(zhuǎn)換成兩兩樣本之間核函數(shù)的計算，大大降低了計算復(fù)雜度和計算量。

1.3 二次互信息及其估計

互信息(mutual information，MI)與二次互信息(quadratic mutual information，QMI)都能夠衡量隨機(jī)變量之間的非線性關(guān)系，但QMI的計算復(fù)雜度小于MI且具有更強(qiáng)的噪聲魯棒性，更適合處理復(fù)雜的工業(yè)過程數(shù)據(jù)。對于兩個隨機(jī)變量X與Y，假設(shè)它們的邊緣概率密度函數(shù)分別為fX(x) 與fY(y)， 聯(lián)合概率密度函數(shù)為fXY(x,y)， 則基于歐氏距離的X與Y之間的QMI定義為

QMIE(X,Y)=?fXY2(x,y)dxdy+?fX2(x)fY2(y)dxdy-2?fXY(x,y)fX(x)fY(y)dxdy

(5)

可以看出，QMIE(X,Y)≥0恒成立。當(dāng)且僅當(dāng)X與Y相互獨立時，QMIE(X,Y)=0。

接下來對QMI進(jìn)行估計。由式(5)，將QMIE(X,Y) 改寫成以下形式

QMIE(X,Y)=VJS+VMS-2VCS

(6)

其中

VJS=?fXY2(x,y)dxdy
VMS=?fX2(x)fY2(y)dxdy
VCS=?fXY(x,y)fX(x)fY(y)dxdy

(7)

選用式(4)中的高斯核函數(shù)，使用Parzen窗估計對VJS、VMS與VCS分別進(jìn)行計算

(8)

2 核二次互信息回歸模型(KQMIR)

2.1 基于核二次互信息(Kernel QMI)的特征提取方法

圖2為基于再生核希爾伯特空間(reproducing kernel hilbert space，RKHS)對輸入X進(jìn)行特征提取的示意圖。首先，將過程數(shù)據(jù)X(輸入)進(jìn)行非線性變換，通過核映射獲取高維特征Φ。之后基于RKHS空間，對高維數(shù)據(jù)進(jìn)行線性變換W提取過程數(shù)據(jù)的低維特征T。低維特征T即為所要求取的過程特征?；谶^程特征T與質(zhì)量數(shù)據(jù)(輸出)Y之間的QMI及過程特征T的Renyi二次熵定義目標(biāo)函數(shù)J，使用梯度下降法求取線性變換W的最優(yōu)值，進(jìn)而求取過程特征T。最后建立過程特征T與質(zhì)量數(shù)據(jù)Y之間的回歸模型進(jìn)行質(zhì)量預(yù)測。

圖2 基于核二次互信息的過程數(shù)據(jù)特征提取

(9)

則不同的核參數(shù)σ將輸入X映射到不同的無窮維的特征空間中，即Q→∞。

接下來，對核矩陣K實施線性變換，將其映射到低維特征空間。首先尋找出一個變換權(quán)重向量w，使得變換后的低維特征t與輸出Y之間的二次互信息最大，同時盡可能保留原始輸入X的自信息。令t=Kw， 其中，w為N×1 權(quán)重向量。t=(t1,…,tN)T為線性變換后的投影值，也是N×1維向量。若共計提取r維新特征，則權(quán)重向量w變?yōu)镹×r維權(quán)重矩陣W，W中的每一列為一個維N×1權(quán)重向量。t變?yōu)镹×r維矩陣T，T的第r列為K在W的第r列wr作用下的得分向量，T為K在低維空間中的新特征，即T=ΦΦTW。 特征提取的目標(biāo)函數(shù)定義如下

J=-HR2(T)·QMIE(T,Y)

(10)

可以看出，通過求取目標(biāo)函數(shù)的最小值，可以使得變換后的過程特征T與輸出Y之間的QMI最大，保證原始輸入數(shù)據(jù)能夠最大程度的包含輸出數(shù)據(jù)所承載的信息量。同時，新提取的T的Renyi二次熵能夠最大程度保留X的自信息。

2.2 模型學(xué)習(xí)過程

(11)

其中

(12)

(13)

由T=KW， 則

(14)

(15)

(16)

(17)

(18)

(19)

(20)

式(9)中核參數(shù)σ的取值可使用輸入樣本X的樣本協(xié)方差矩陣對角線元素σx確定。使用Silverman’s經(jīng)驗準(zhǔn)則計算σt，即

(21)

其中，r為特征空間的維數(shù)。σy可以直接由輸出樣本的方差確定。

2.3 回歸預(yù)測模型

經(jīng)過以上計算，由梯度下降法求得W*，此時T*=ΦΦTW*。 建立T與Y之間的最小二乘回歸模型為

(22)

將T*=ΦΦTW*代入式(22)，則X與Y之間的最終回歸模型為

(23)

θ為回歸系數(shù)

θ=ΦTW*(W*TK2W*)-1W*TKY

(24)

將以上模型稱為二次互信息回歸(kernel quadratic mutual information regression，KQMIR)模型。

3 基于KQMIR模型的質(zhì)量預(yù)測

3.1 數(shù)據(jù)預(yù)處理

發(fā)酵過程是間歇過程，在一個批次的生產(chǎn)結(jié)束后才能獲取產(chǎn)品。因此，發(fā)酵過程的可測量過程數(shù)據(jù)是三維信息，包括時間、變量與批次。實際生產(chǎn)中，會記錄當(dāng)前批次的產(chǎn)品質(zhì)量變量信息，比如菌體濃度等。通常將三維矩陣按照建模的需要按不同方式展開成二維矩陣[9]。本文選用批次展開方式，得到二維矩陣X(I×KJx)， 如圖3所示。其中，I為批次數(shù)，Jx為過程變量個數(shù)，K為總采樣時刻。同時，每個批次結(jié)束后，會產(chǎn)生相應(yīng)的產(chǎn)品質(zhì)量變量，因此，質(zhì)量數(shù)據(jù)陣可表示為Y(I×L)，L為質(zhì)量變量的個數(shù)。

圖3 間歇過程三維數(shù)據(jù)的批次展開方式

3.2 離線建模

離線建模步驟如下：

步驟2 對X與Y沿批次方向進(jìn)行中心化與方差歸一化，求取核參數(shù)σ及σx，σt，σy；

步驟3 使用核函數(shù)求取核矩陣Ko，并對其進(jìn)行中心化及方差歸一化。中心化處理如下

(25)

其中，IN是N維單位矩陣，1N是元素全為1的N×1維向量。方差歸一化處理后得到最終K如下

(26)

步驟4 初始化權(quán)重矩陣W(N×r)。 使用梯度下降法求取目標(biāo)函數(shù)J的最小值，計算出最優(yōu)投影W*。梯度下降法的迭代過程如下：

(1)使用式(3)至式(8)計算J；

(3)使用式(11)更新權(quán)重矩陣，保證W的列向量之間相互正交；

(4)使用T=ΦΦTW計算特征T；

(5)重復(fù)以上過程，直到達(dá)到終止條件。

3.3 在線質(zhì)量預(yù)測

在線質(zhì)量預(yù)測步驟如下：

步驟1 對當(dāng)前第k時刻樣本xk(Jx×1) 進(jìn)行數(shù)據(jù)填充。下一時刻至第K個采樣時刻的樣本值均使用xk填充，得到新樣本xnew(KJx×1)=xnew(M×1)。 依據(jù)X的均值和標(biāo)準(zhǔn)差對xnew進(jìn)行中心化與方差歸一化處理；

步驟2 計算xnew的核向量Knew

Knew=φT(xnew)·ΦT=
[φT(xnew·φ(x1)…φT(xnew)·φ(xi)…φT(xnew)·φ(xN))]

(27)

對其中心化如下

(28)

方差歸一化后得到最終Knew，計算如下

(29)

步驟3 依據(jù)式(23)～式(24)使用KQMIR模型進(jìn)行質(zhì)量預(yù)測，當(dāng)前時刻對質(zhì)量的預(yù)測值為

(30)

步驟5 使用以上步驟1至步驟4順序計算未來每個采樣時刻的質(zhì)量預(yù)測值，直到當(dāng)前批次結(jié)束。

最后，基于核二次互信息回歸模型的質(zhì)量預(yù)測的總體流程如圖4所示。

圖4 基于KQMIR的在線質(zhì)量預(yù)測

本文采用批次RMSE(簡稱RMSE)與采樣時刻RMSE(簡稱RMSE(k))評價模型預(yù)測效果。批次RMSE定義如下

(31)

(32)

4 實驗及結(jié)果分析

4.1 青霉素發(fā)酵仿真實驗

Pensim是國際上較為廣泛應(yīng)用的青霉素發(fā)酵過程仿真平臺。使用Pensim平臺產(chǎn)生48批次的正常生產(chǎn)數(shù)據(jù)，任選40批次用于建模，8批次進(jìn)行測試。發(fā)酵周期為400 h，采樣間隔為1 h。共計選擇10個過程變量(通風(fēng)速率、攪拌功率、冷水流加速率、反應(yīng)器體積、產(chǎn)熱、溶氧濃度、溫度、補(bǔ)料溫度、pH、排氣CO2濃度)與一個質(zhì)量變量(青霉素濃度, 單位g/L)。

特征維數(shù)r不同取值下，KQMIR對8批次測試數(shù)據(jù)的批次RMSE的均值見表1。當(dāng)r=4時，KQMIR預(yù)測誤差最小。圖5為不同r取值下KQMIR的采樣時刻RMSE(k)，200 h-400 h區(qū)間內(nèi)，當(dāng)r取2或4時，RMSE(k)取值較小，當(dāng)r=1時，RMSE(k)取值較大，與表1結(jié)果一致。綜上，在r=4時預(yù)測效果最優(yōu)。

表1 Pensim中KQMIR模型的批次RMSE均值

圖5 KQMIR對所有測試批次的RMSE(k)

為驗證方法的有效性，將本文提出的KQMIR模型與核偏最小二乘法(kernel partial least squares，KPLS)及支持向量回歸法(support vector regression，SVR)進(jìn)行對比實驗。圖6為某個測試批次下3種方法對批次最終青霉素濃度的預(yù)測結(jié)果。在前200 h，KQMIR與KPLS的預(yù)測效果明顯優(yōu)于SVR，且KQMIR預(yù)測效果略優(yōu)于KPLS。在批次結(jié)束時，KQMIR方法預(yù)測結(jié)果最接近測量值。表2列出3種方法對所有8個測試批次的RMSE的均值及標(biāo)準(zhǔn)差。由RMSE的均值可以看出KQMIR對青霉素濃度的平均預(yù)測誤差小于KPLS及SVR。從RMSE的標(biāo)準(zhǔn)差可以看出，SVR對不同測試批次的預(yù)測效果最穩(wěn)定，KQMIR的預(yù)測穩(wěn)定性略優(yōu)于KPLS。圖7為3種方法對所有測試批次的RMSE(k)。在生產(chǎn)開始及210 h前的絕大多數(shù)采樣時刻，SVR的RMSE(k)值波動幅度較大，且遠(yuǎn)大于KQMIR及KPLS。210 h后，SVR的RMSE(k)最小，KQMIR的RMSE(k)值小于KPLS。

圖7 3種方法對所有測試批次的RMSE(k)

表2 3種方法的批次RMSE均值及標(biāo)準(zhǔn)差(Pensim)

圖6 3種方法對某個測試批次的青霉素濃度預(yù)測結(jié)果

綜上，KQMIR對青霉素濃度的平均預(yù)測誤差最小，預(yù)測效果最好。與在高維特征中采用相關(guān)性作為準(zhǔn)則提取特征的KPLS方法相比，KQMIR在高維特征空間中采用二次互信息作為準(zhǔn)則自適應(yīng)提取特征。由于二次互信息是更高階的統(tǒng)計量，因此能夠更好地提取原始數(shù)據(jù)所包含的信息。

4.2 大腸桿菌發(fā)酵生產(chǎn)數(shù)據(jù)實驗

實際生產(chǎn)數(shù)據(jù)驗證更突顯本方法的意義和效果，本文以北京亦莊某制藥公司使用大腸桿菌制備抗腫瘤藥用蛋白白介素-2(IL-2)的生產(chǎn)數(shù)據(jù)為驗證對象。實際生產(chǎn)中，產(chǎn)物IL-2的質(zhì)量往往難以在線預(yù)測，需要實驗室分析。本文將對IL-2質(zhì)量進(jìn)行在線預(yù)測。

實驗中，共計選擇8個可測量的生產(chǎn)過程變量(pH值、罐壓、溶氧、攪拌速率、溫度、通風(fēng)速率、補(bǔ)氮、補(bǔ)碳)與一個質(zhì)量變量(IL-2質(zhì)量,單位kg)。發(fā)酵時長約為19 h-20 h，采樣間隔為0.5 h，共計39個采樣時刻。共收集到28批次正常生產(chǎn)數(shù)據(jù)，隨機(jī)選取其中的20批次進(jìn)行離線回歸建模，剩余8批次用于測試。實驗開發(fā)語言為Matlab，開發(fā)平臺為Matlab R2014a。

特征維數(shù)r取不同值時，KQMIR對8批次測試數(shù)據(jù)的預(yù)測RMSE的均值見表3，可以看出當(dāng)r=2時，KQMIR取得最佳預(yù)測效果，具有最小的RMSE。此外，不同r取值下KQMIR的RMSE(k)值如圖8所示，可以看出，RMSE(k)值大致在1-1.5區(qū)間內(nèi)波動，在采樣時刻22之后直至批次結(jié)束這段區(qū)間內(nèi)，r=2時都具有最小的RMSE(k)。

表3 KQMIR模型的批次RMSE均值

圖8 KQMIR對所有測試批次的RMSE(k)

圖9為某個測試批次下3種方法對IL-2質(zhì)量的預(yù)測結(jié)果對比?？梢钥闯?，KQMIR的預(yù)測效果最好、最接近測量值。在采樣時刻25至批次結(jié)束期間，KQMIR的預(yù)測結(jié)果好于KPLS，明顯優(yōu)于SVR。

圖9 3種方法對某個測試批次IL-2質(zhì)量的預(yù)測結(jié)果

3種方法對所有測試批次的RMSE的均值及標(biāo)準(zhǔn)差見表4，由RMSE均值可以看出KQMIR對IL-2質(zhì)量的預(yù)測精度優(yōu)于KPLS，且明顯優(yōu)于SVR。另一方面，從RMSE的標(biāo)準(zhǔn)差可以看出，SVR對所有測試批次的預(yù)測效果最穩(wěn)定，KQMIR的預(yù)測穩(wěn)定性要略優(yōu)于KPLS。此外，3種方法對所有測試批次的RMSE(k)值如圖10所示，可以看出，SVR的RMSE(k)值明顯大于KQMIR及KPLS。采樣時刻15以后的所有采樣時刻上，KQMIR的RMSE(k)值均小于KPLS。綜上，本文提出的KQMIR具有最優(yōu)的預(yù)測效果。

表4 3種方法的批次RMSE均值及標(biāo)準(zhǔn)差

圖10 3種方法對所有測試批次的RMSE(k)

5 結(jié)束語

為解決發(fā)酵過程數(shù)據(jù)的強(qiáng)非線性、質(zhì)量變量難以在線預(yù)測的問題，提出一種基于二次互信息回歸模型的質(zhì)量預(yù)測方法。該方法具有以下優(yōu)點：①將原始輸入空間通過核映射到高維空間，加強(qiáng)對非線性數(shù)據(jù)的處理能力；②所提取的過程特征能夠最好地描述過程數(shù)據(jù)和質(zhì)量數(shù)據(jù)之間互相包含的信息，盡可能使得過程變量和質(zhì)量變量之間的二次互信息最大化；③提取的特征可以最大化地保留過程數(shù)據(jù)的自信息量，盡可能使得過程變量自身的Renyi二次熵最大。該方法與經(jīng)典的非線性回歸方法KPLS及SVR進(jìn)行對比實驗。青霉素發(fā)酵仿真實驗中對青霉素濃度的預(yù)測及大腸桿菌發(fā)酵過程中對IL-2質(zhì)量進(jìn)行預(yù)測的實驗，驗證了KQMIR的有效性，并且KQMIR具有優(yōu)于KPLS及SVR的預(yù)測效果。KQMIR的質(zhì)量預(yù)測方法可以推廣到其它間歇過程中。目前，二次互信息準(zhǔn)則還主要應(yīng)用在圖像分類、圖像識別等領(lǐng)域，在過程建模中很少應(yīng)用。本文將其引入過程領(lǐng)域，旨在探索QMI的相關(guān)理論在過程數(shù)據(jù)處理中的應(yīng)用。