蔡子君，謝 莉，楊慧中

(江南大學(xué) 輕工過程先進(jìn)控制教育部重點(diǎn)實(shí)驗室，江蘇 無錫 214122)

0 引言

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

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

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

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

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

4)測量難度大。許多關(guān)鍵變量無法在線測量，需要進(jìn)行離線分析，例如青霉素濃度、氮濃度、濃稠度等，導(dǎo)致關(guān)鍵變量的數(shù)據(jù)稀缺。

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

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

1 青霉素發(fā)酵階段劃分

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

圖2 工業(yè)環(huán)境下冷水流加速率圖

考慮實(shí)際工業(yè)過程的操作環(huán)境，本文將使用葡萄糖流加速率作為調(diào)度變量。葡萄糖是青霉素發(fā)酵中底物的一種，其流加速率雖然在發(fā)酵過程中經(jīng)過人為調(diào)整，但調(diào)整頻率相對較低并且整體曲線能夠體現(xiàn)發(fā)酵過程的階段特性[10]。

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

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

2 模型結(jié)構(gòu)：FIR融合模型

考慮到在工業(yè)現(xiàn)場青霉素濃度為慢速采樣[13]，無法用于模型輸入，所以本文采用FIR模型作為局部模型，然后通過權(quán)重函數(shù)將3個局部模型融合為全局模型來描述青霉素發(fā)酵的動態(tài)特性。局部FIR模型結(jié)構(gòu)如下：

(1)

(2)

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

采用高斯分布作為權(quán)重函數(shù)[14]：

(3)

(4)

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

(5)

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

3 基于VB算法的模型參數(shù)辨識

3.1 VB算法回顧

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

(6)

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

(7)

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

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

(8)

將對數(shù)邊緣函數(shù)與下界函數(shù)做差，可得：

KL(q|p)

(9)

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

3.2 基于VB的辨識算法推導(dǎo)

3.2.1 模型參數(shù)的先驗分布

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

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

將上述先驗分布用一個聯(lián)合先驗分布表示：

(10)

3.2.2 VB算法：E步

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

(11)

求解如下優(yōu)化問題：

計算關(guān)于q(I)的拉格朗日函數(shù)的導(dǎo)數(shù)：

逐項求導(dǎo)后可以得到：

(12)

令：

得到：

(13)

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

(14)

由此可以將ZI表達(dá)為：

繼而定義：

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

(15)

為簡化表達(dá)，令：

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

式(15)可簡化為：

(16)

(17)

其中:

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

3.2.3 VB算法：M步

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

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

(18)

其中：

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

1)q(θi)部分：計算關(guān)于q(θi)的拉格朗日函數(shù)的一階導(dǎo)數(shù)

可得:

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

即:

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

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

其中：

(19)

2)q(ηi)部分：類似地關(guān)于q(ηi)對式(18)的求導(dǎo):

可得:

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

其中:

整理后得到:

(20)

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

3)q(σi)部分：關(guān)于q(σi)對式(18)的三部分求導(dǎo):

可得:

(21)

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

根據(jù)伽馬分布的相關(guān)知識，可得:

(22)

(23)

最后通過非線性數(shù)值優(yōu)化的方法求得局部模型寬度Oi的點(diǎn)估計，優(yōu)化目標(biāo)函數(shù)如下:

(24)

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

3.2.4 VB辨識算法計算步驟

基于VB的參數(shù)辨識算法計算步驟總結(jié)如下。

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

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

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

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

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

4 仿真驗證

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

圖4 預(yù)處理前后排放氣體中CO2濃度圖

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

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

圖5 青霉素濃度擬合曲線

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

圖6 對比實(shí)驗結(jié)果圖

可以看出，通過對青霉素發(fā)酵準(zhǔn)確的階段劃分并充分考慮工業(yè)環(huán)境中的變化因素，本文所建立的模型能夠更好地預(yù)測實(shí)際工業(yè)環(huán)境中青霉素發(fā)酵過程的產(chǎn)物濃度，對青霉素生產(chǎn)過程的控制與優(yōu)化具有一定的指示作用。

5 結(jié)束語

本文采用變分貝葉斯算法建立了基于不同工作點(diǎn)的青霉素發(fā)酵過程各階段子模型，采用由調(diào)度變量計算得到的標(biāo)準(zhǔn)化權(quán)值將子模型進(jìn)行融合，變分貝葉斯算法通過估計參數(shù)的后驗分布，能夠?qū)l(fā)酵中過程變量的不確定利用均值、方差等統(tǒng)計特性解析地表達(dá)出來，在環(huán)境復(fù)雜的工業(yè)青霉素發(fā)酵過程中表現(xiàn)出優(yōu)異的性能。仿真實(shí)驗表明本文方法能夠精確建立青霉素發(fā)酵過程產(chǎn)物濃度模型，在復(fù)雜的工業(yè)環(huán)境中準(zhǔn)確地預(yù)測青霉素發(fā)酵過程。