基于變分貝葉斯算法的青霉素發酵過程建模

蔡子君，謝 莉，楊慧中

(江南大學 輕工過程先進控制教育部重點實驗室，江蘇 無錫 214122)

0 引言

抗生素是一類生物產生的具有抑制某些細胞生長的次級代謝產物[1]，其對于某些病原微生物的抑制和滅殺作用使其成為一種重要的藥物。作為一種典型的抗生素，青霉素發酵過程具有以下特點：

1)時變性。發酵過程中青霉素的濃度取決于生成青霉素的菌絲濃度，而影響菌絲生長的因素眾多，例如發酵罐中糖濃度、溶解氧濃度、pH值、溫度等，并且這些變量都會隨著時間改變，導致青霉素的發酵過程呈現很強的時變性。

2)非線性。青霉素是一種次級代謝產物，生成青霉素所需的眾多基質、前體、以及副產物之間會發生復雜的反應，形成了一個多輸入多輸出的非線性系統，加大了青霉素的濃度預測難度。

3)階段性。青霉素發酵通常經歷4個階段，分別為準備期、對數生長期、平穩期和消亡期[2]，如圖1所示，準備期中菌絲慢慢生成，青霉素主要在對數生長期生成，經過平穩期后菌絲逐漸水解，進入消亡期。

圖1 青霉素濃度曲線階段圖

4)測量難度大。許多關鍵變量無法在線測量，需要進行離線分析，例如青霉素濃度、氮濃度、濃稠度等，導致關鍵變量的數據稀缺。

上述特點使得對青霉素發酵過程的預測和控制較為困難。過去提出的各種青霉素發酵過程的預測模型大致可以分為機理模型和數據驅動模型兩類。最著名的機理模型是由Birol于2002年改進的非結構模型[3]。由于青霉素化學結構復雜，發酵過程涉及的反應眾多，大部分機理模型不對內部結構討論[4]，而是對青霉素與菌絲的不同部分之間的動態關系進行描述。這種非結構模型能夠一定程度上描述反應過程并在特定參數環境下模擬發酵過程，但是在面對實際發酵過程多變的環境時就顯露出適應性差的缺點。近年來，許多基于數據驅動的青霉素發酵模型被提出，利用諸如神經網絡、支持向量機等方法模擬青霉素發酵過程[5]。相對機理模型，數據驅動模型具有以下兩點優勢：1)不需要復雜的發酵機理知識，降低了建模難度；2)在發酵環境改變時不需要大量實驗重新確定模型參數，降低了模型的維護成本。但現有的數據驅動模型大多利用Pensim仿真平臺產生的數據，只能描述出在實驗室規模下青霉素發酵的大致過程。

本文針對實際工業中的青霉素發酵過程，基于變分貝斯算法建立FIR融合模型。在模型選擇方面，由于青霉素發酵過程具有明顯的階段性[6]，本文采用多模型融合的思路，選取能夠反應階段特征的關鍵過程變量作為調度變量進行階段劃分，確定幾個典型工況點以計算各局部模型的權重。工業青霉素發酵過程中生產環境復雜性使得過程存在著不確定性，這對發酵過程的辨識造成了很大的困難[7]。不確定性在工業過程辨識中的處理方法可分為兩種[8]，一種是通過不同模型的變化體現系統的不確定性；而本文采用的是第二種，即在模型結構已知的條件下，將系統不確定性通過模型參數的不確定性來體現。本文采用變分貝葉斯算法作為辨識算法，它是期望最大化(EM)算法在貝葉斯方法上的一種推廣，相較于EM算法的點估計，變分貝葉斯算法可以估計參數和隱變量的整個后驗概率分布，能夠更好地描述青霉素發酵過程的不確定性。

1 青霉素發酵階段劃分

調度變量是指可以反應系統工作狀態與階段特征并能夠人為調控的過程變量[9]。在Pensim仿真平臺的環境下，冷水流加速率能比較好的反映發酵的階段特征，所以在以往的青霉素發酵模型中經常被作為模型的調度變量。但是，在實際工業過程中，操作人員會頻繁地改變冷水流加速率以保持穩定的發酵罐溫度，導致實際過程中冷水流加速率(圖2)不宜作為過程的調度變量。

圖2 工業環境下冷水流加速率圖

考慮實際工業過程的操作環境，本文將使用葡萄糖流加速率作為調度變量。葡萄糖是青霉素發酵中底物的一種，其流加速率雖然在發酵過程中經過人為調整，但調整頻率相對較低并且整體曲線能夠體現發酵過程的階段特性[10]。

為劃分發酵過程的各個階段，需要進一步對調度變量進行聚類。本文采用模糊C均值(Fuzzy C-Means，FCM)聚類方法[11]對葡萄糖流加速率進行聚類，并將各聚類中心作為發酵過程的典型工作點。FCM是一種基于目標函數的聚類算法，首先對聚類中心設置初值，然后通過最小化數據點與各聚類中心的距離和模糊隸屬度的加權和為目標，不斷修正聚類中心和分類矩陣直到符合終止準則，將具有類似特征的數據聚為一類。

在青霉素發酵過程的準備期中，菌體快速生長消耗氧氣，當溶解氧水平下降到一定程度時，菌體會產生中間代謝物并開始生成青霉素[12]，所以本文假定在準備期中青霉素濃度為零，而對數生長期、平穩期和消亡期這3個階段分別對應3個不同的工作點。因此，在采用FCM算法對實際工業青霉素發酵過程中的葡萄糖流加速率進行聚類時，將聚類中心個數設置為3，相應的聚類結果如圖3所示，其中三角形表示計算得到的聚類中心。

2 模型結構：FIR融合模型

考慮到在工業現場青霉素濃度為慢速采樣[13]，無法用于模型輸入，所以本文采用FIR模型作為局部模型，然后通過權重函數將3個局部模型融合為全局模型來描述青霉素發酵的動態特性。局部FIR模型結構如下：

(1)

(2)

其中：k=1,2…N表示發酵過程的各時刻,輸出變量yk為k時刻發酵罐中青霉素濃度，Ik=1, 2, 3表示k時刻變量隸屬的局部模型，模型選擇的輸入變量個數為m，輸入階數設為na，ek是均值為0方差為σi-1(未知的待辨識參數)的高斯分布噪聲，下標i=Ik表示各個局部模型。

采用高斯分布作為權重函數[14]：

(3)

(4)

若各個局部模型的有效寬度Oi已知，根據式(3)～(4)計算出各子模型的權重后，容易得到系統全局模型的預測輸出：

(5)

然而，局部模型的有效寬度Oi未知，需要與各子模型的參數θi以及噪聲方差σi-1同時辨識，增加了系統辨識的難度。此外，考慮到實際過程的不確定性，本文通過引入參數不確定性，在VB算法框架下推導相應的辨識算法。

3 基于VB算法的模型參數辨識

3.1 VB算法回顧

令模型的參數集為Θ，隱變量集為Cmis，觀測數據集為Cobs。對于存在未知參數的模型，邊緣似然函數可以由下式計算：

(6)

而式(6)中含有難以計算的高維積分，VB算法通過構造聯合分布q(Cmis,Θ)來近似計算后驗分布p(Cmis,Θ)，運用Jensen不等式[15]得到：

(7)

假定聯合分布q(Cmis,Θ)是可分解的[16]，得到對數邊緣函數的下界函數：

F[q(Cmis)q(Θ)]=

(8)

將對數邊緣函數與下界函數做差，可得：

KL(q|p)

(9)

與期望最大化算法類似，變分貝葉斯算法不斷迭代地更新隱變量和模型參數的后驗分布，直至算法收斂，得到真實后驗分布p(Cmis,Θ|Cobs)的近似分布q(Cmis,Θ)[18]。

3.2 基于VB的辨識算法推導

3.2.1 模型參數的先驗分布

p(θi|ηi) =N(0,ηi-1Dna ×na)

其中:D表示單位矩陣，g表示伽馬分布,a0,b0,c0,d0為伽馬分布的超參數,Γ表示伽馬函數，其表達式為:

將上述先驗分布用一個聯合先驗分布表示：

(10)

3.2.2 VB算法：E步

在E步驟中，固定參數集，關于隱變量對下界函數求極值，得到隱變量的更新后驗分布q(I)。下界函數F[q(I),q(Θ)]可表示為:

(11)

求解如下優化問題：

計算關于q(I)的拉格朗日函數的導數：

逐項求導后可以得到：

(12)

令：

得到：

(13)

其中:<·>q(Θ)代表對q(Θ)求期望。

(14)

由此可以將ZI表達為：

繼而定義：

其中:q(Ik=i)可以表示為:

(15)

為簡化表達，令：

I1:N|Θ,O)>q(Θ)}

式(15)可簡化為：

(16)

(17)

其中:

將Aki代入式(16)即可得到q(Ik=i)。

3.2.3 VB算法：M步

在M步驟中，固定隱變量，關于參數集對下界函數(11)求極值，得到參數集的更新式q(Θ)。下界函數可以表示為：

F[q(I),q(Θ)]=

(18)

其中：

下面利用拉格朗日乘子法，依次關于q(θi),q(ηi)和q(σi)最大化下界函數。

1)q(θi)部分：計算關于q(θi)的拉格朗日函數的一階導數

可得:

-lnq(θi)+lnp(θi|<η1>q(η))+Cθ-

即:

q(θi)∝p(θi|<ηi>q(η))×

因此，q(θi)服從高斯分布，即:

其中：

(19)

2)q(ηi)部分：類似地關于q(ηi)對式(18)的求導:

可得:

-lnq(ηi)+(a0+1)lnp(ηi)-

其中:

整理后得到:

(20)

其中:Cη是與ηi無關的常數。由式(20)可知q(ηi)服從伽馬分布，即q(ηi)=g(ai,bi)，且:

3)q(σi)部分：關于q(σi)對式(18)的三部分求導:

可得:

(21)

其中:Cσ是與σi無關的常數。根據式(21)可知q(σi)服從伽馬分布，即q(σi)=g(ci,di)，且:

根據伽馬分布的相關知識，可得:

(22)

(23)

最后通過非線性數值優化的方法求得局部模型寬度Oi的點估計，優化目標函數如下:

(24)

其中:q(Ik=i)由式(16)計算得到:

3.2.4 VB辨識算法計算步驟

基于VB的參數辨識算法計算步驟總結如下。

1)根據式(10)給模型的未知參數Θ設置合適的先驗分布p(Θ)，初始化未知參數Θ、O以及先驗分布中的超參數a0,b0,c0,d0。

2)E步：關于隱變量最大化下界函數，根據式(17)計算Aki，并由式(16)得到隱變量Ik后驗分布的更新式q(Ik)。

3)M步：固定隱變量，分別關于q(θi)、q(ηi)、q(σi)最大化下界函數，根據式(19)、式(22)以及式(23)得到各參數的期望E(θ)、E(η)以及E(σ)，并根據式(24)計算局部模型有效寬度Oi的點估計。

4)根據以上兩步得到的模型參數以及隱變量計算下界函數F：

5)將獲得的估計值代入2)，不斷迭代計算2)～4)，直至下界函數收斂。

4 仿真驗證

本文利用來自文獻[13]的10個100 000 L的青霉素流加發酵罐產生的數據訓練并測試模型以驗證模型對實際工業環境的適應性。實際工業過程中采集到的數據夾雜著大量噪聲，以發酵罐排放氣體中的CO2濃度(圖4)為例。

圖4 預處理前后排放氣體中CO2濃度圖

為減少數據中的噪聲對模型穩定性的影響，首先對輸入數據中的異常值進行處理，由于發酵時間較長導致數據前后浮動較大，故本文首先根據發酵過程的階段特性分階段利用3σ準則剔除了輸入數據中的異常值，然后對輸入作濾波處理。

輸入變量選擇方面，在充分考慮反應機理并計算各輸入變量與青霉素濃度相關系數后，選取5個過程變量作為輸入變量，分別為：排放氣體中CO2百分比(%)，排放氣體中O2百分比(%)，pH值，C的生成率(g/min)，發酵罐中物質總質量(kg)。利用本文第1節對葡萄糖流加速率進行聚類得到的典型工作點和經過預處理的輸入數據，應用本文所述方法辨識得到子模型的參數，融合后得到青霉素濃度擬合曲線，如圖5所示，其中三角形表示實際測量得到的數據。模型的質量通過相關誤差進行測量，其計算公式如下：

圖5 青霉素濃度擬合曲線

選取正常發酵情況下的另一批次數據對得到的模型進行測試，并與文獻[9]中基于EM算法的青霉素發酵建模方法進行對比，預測結果如圖6所示，計算得到本文模型與文獻[9]的相關誤差分別為0.23%和0.75%。

圖6 對比實驗結果圖

可以看出，通過對青霉素發酵準確的階段劃分并充分考慮工業環境中的變化因素，本文所建立的模型能夠更好地預測實際工業環境中青霉素發酵過程的產物濃度，對青霉素生產過程的控制與優化具有一定的指示作用。

5 結束語

本文采用變分貝葉斯算法建立了基于不同工作點的青霉素發酵過程各階段子模型，采用由調度變量計算得到的標準化權值將子模型進行融合，變分貝葉斯算法通過估計參數的后驗分布，能夠將發酵中過程變量的不確定利用均值、方差等統計特性解析地表達出來，在環境復雜的工業青霉素發酵過程中表現出優異的性能。仿真實驗表明本文方法能夠精確建立青霉素發酵過程產物濃度模型，在復雜的工業環境中準確地預測青霉素發酵過程。