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

李 征，王 普，高學金+，齊詠生，高慧慧

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

0 引 言

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

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

1 相關理論

1.1 發(fā)酵過程工藝

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

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

1.2 Renyi二次熵及其估計

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

(1)

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

(2)

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

(3)

(4)

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

1.3 二次互信息及其估計

互信息(mutual information，MI)與二次互信息(quadratic mutual information，QMI)都能夠衡量隨機變量之間的非線性關系，但QMI的計算復雜度小于MI且具有更強的噪聲魯棒性，更適合處理復雜的工業(yè)過程數據。對于兩個隨機變量X與Y，假設它們的邊緣概率密度函數分別為fX(x) 與fY(y)， 聯合概率密度函數為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恒成立。當且僅當X與Y相互獨立時，QMIE(X,Y)=0。

接下來對QMI進行估計。由式(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)中的高斯核函數，使用Parzen窗估計對VJS、VMS與VCS分別進行計算

(8)

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

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

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

圖2 基于核二次互信息的過程數據特征提取

(9)

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

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

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

(10)

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

2.2 模型學習過程

(11)

其中

(12)

(13)

由T=KW， 則

(14)

(15)

(16)

(17)

(18)

(19)

(20)

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

(21)

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

2.3 回歸預測模型

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

(22)

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

(23)

θ為回歸系數

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

(24)

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

3 基于KQMIR模型的質量預測

3.1 數據預處理

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

圖3 間歇過程三維數據的批次展開方式

3.2 離線建模

離線建模步驟如下：

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

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

(25)

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

(26)

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

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

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

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

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

3.3 在線質量預測

在線質量預測步驟如下：

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

步驟2 計算xnew的核向量Knew

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

(27)

對其中心化如下

(28)

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

(29)

步驟3 依據式(23)～式(24)使用KQMIR模型進行質量預測，當前時刻對質量的預測值為

(30)

步驟5 使用以上步驟1至步驟4順序計算未來每個采樣時刻的質量預測值，直到當前批次結束。

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

圖4 基于KQMIR的在線質量預測

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

(31)

(32)

4 實驗及結果分析

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

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

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

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

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

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

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

表2 3種方法的批次RMSE均值及標準差(Pensim)

圖6 3種方法對某個測試批次的青霉素濃度預測結果

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

4.2 大腸桿菌發(fā)酵生產數據實驗

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

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

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

表3 KQMIR模型的批次RMSE均值

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

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

圖9 3種方法對某個測試批次IL-2質量的預測結果

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

表4 3種方法的批次RMSE均值及標準差

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

5 結束語

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