劉重陽(yáng)， 尹 蕾， 馮恩民， 修志龍

（1.大連理工大學(xué)數(shù)學(xué)科學(xué)學(xué)院，遼寧大連 116024；2.山東工商學(xué)院數(shù)學(xué)與信息科學(xué)學(xué)院，山東煙臺(tái) 264005；3.大連理工大學(xué)生命科學(xué)與技術(shù)學(xué)院，遼寧大連 116024）

0 引 言

1，3－丙二醇（1，3－PD）是一種重要的化工原料，可用來(lái)合成許多具有優(yōu)良特性的聚合物，如聚酯和聚氨酯［1］.傳統(tǒng)的1，3－丙二醇化學(xué)合成法生產(chǎn)需要高溫、高壓及貴重催化劑才能實(shí)現(xiàn)，此法分離提純困難、成本高，極大地限制了1，3－PD的發(fā)展.隨著現(xiàn)代生物技術(shù)的發(fā)展，人們開(kāi)始嘗試用微生物發(fā)酵法生產(chǎn)1，3－PD.微生物發(fā)酵法具有條件溫和、操作簡(jiǎn)單、副產(chǎn)物少、綠色環(huán)保等優(yōu)點(diǎn)，所以這方面的研究正受到國(guó)內(nèi)外越來(lái)越多的重視.特別地，通過(guò)克雷伯氏桿菌（K.pneumoniae）發(fā)酵生產(chǎn)1，3－PD受到眾多學(xué)者的廣泛關(guān)注［2～5］.利用克雷伯氏桿菌進(jìn)行甘油發(fā)酵是一個(gè)復(fù)雜的生化過(guò)程，因?yàn)槲⑸锏纳L(zhǎng)受到底物、各種產(chǎn)物和中間代謝物的抑制作用［3］，其中3－羥基丙醛（3－HPA）是一種有毒的中間代謝物，它對(duì)生產(chǎn)1，3－丙二醇所必需的兩種重要的酶（甘油脫水酶（GDHt）和1，3丙二醇氧化還原酶（PDOR））的活性具有抑制作用.近年來(lái)，國(guó)內(nèi)外許多研究者利用數(shù)學(xué)模型來(lái)描述甘油發(fā)酵過(guò)程以實(shí)現(xiàn)發(fā)酵過(guò)程的優(yōu)化從而提高1，3－PD的產(chǎn)量.Zeng等［3］提出了關(guān)于底物消耗和產(chǎn)物形成的過(guò)量動(dòng)力學(xué)模型.修志龍等［4］改進(jìn)了過(guò)量動(dòng)力學(xué)模型使其適用于更大范圍的甘油注入濃度并研究了連續(xù)發(fā)酵中出現(xiàn)的多穩(wěn)態(tài)現(xiàn)象.但是，上述模型都沒(méi)有考慮胞內(nèi)1，3－PD、胞內(nèi)甘油以及3－HPA的濃度變化.然而，3－HPA對(duì)于1，3－PD的生成起著至關(guān)重要的作用，所以上述模型不能恰當(dāng)?shù)孛枋龈视桶l(fā)酵過(guò)程.

2008年，Sun等首次提出涉及胞內(nèi)1，3－PD跨膜運(yùn)輸及3－HPA的數(shù)學(xué)模型來(lái)描述甘油發(fā)酵生產(chǎn)1，3－PD過(guò)程［6］.但是該模型中的胞內(nèi)動(dòng)力學(xué)參數(shù)的確定是根據(jù)大腸桿菌和其他菌種類(lèi)推得到的，而沒(méi)有利用實(shí)驗(yàn)數(shù)據(jù)及發(fā)酵過(guò)程的信息.顯然這種方法存在一定的偏差.為了更好地描述甘油發(fā)酵生產(chǎn)1，3－PD這一生化過(guò)程，本文改進(jìn)文獻(xiàn)［6］中的數(shù)學(xué)模型，以計(jì)算值與實(shí)驗(yàn)數(shù)據(jù)之間的平均相對(duì)誤差作為優(yōu)化目標(biāo)，建立確定胞內(nèi)動(dòng)力學(xué)參數(shù)的參數(shù)辨識(shí)模型.參數(shù)辨識(shí)對(duì)于確定未知?jiǎng)恿W(xué)方程，進(jìn)一步解釋和定性地描述發(fā)酵行為是非常重要的.目前，關(guān)于參數(shù)辨識(shí)問(wèn)題的求解已有許多算法，其中包括非線性最小二乘法［7］、單純形法［8、9］、遺傳算法［10］和差異演化算法［11］等.本文將粒子群算法和Hooke－Jeeves算法相結(jié)合構(gòu)造改進(jìn)的粒子群算法求解參數(shù)辨識(shí)模型.

1 連續(xù)發(fā)酵動(dòng)力學(xué)系統(tǒng)

在文獻(xiàn)［6］的基礎(chǔ)上，本文建立如下描述甘油連續(xù)發(fā)酵過(guò)程的動(dòng)力學(xué)系統(tǒng)：

式中：t∈［0，T］，T∈（0，＋∞）表示連續(xù)發(fā)酵達(dá)到穩(wěn)態(tài)的時(shí)刻；x1（t），x2（t），…，x8（t）分別表示生物量、胞外甘油、胞外1，3－PD、乙醇、乙酸、胞內(nèi)甘油、3－HPA和胞內(nèi)1，3－PD在t時(shí)刻的濃度；x i（0）（i∈I8）為發(fā)酵的初始濃度；cs0為注入甘油濃度；D表示稀釋速率；Km1、Km2分別表示GDHt和PDOR的米氏常數(shù)，且取值分別為0.53、0.14 mmol／L.基于文獻(xiàn)［4］和［6］，細(xì)胞的比生長(zhǎng)速率μ、比消耗速率p1、比生成速率p2和p3，及u1、u2分別由下式給出：

對(duì)于系統(tǒng)（8），本文以文獻(xiàn)［6］中的參數(shù)值作為辨識(shí)的初始值k0（文獻(xiàn)［6］中無(wú)k1，根據(jù)生物意義此處取值為100；為了保證系統(tǒng)達(dá)到穩(wěn)態(tài)，參數(shù)k3、k5、k10、k11、k15的值作了相應(yīng)的修改），定義參數(shù)向量的允許范圍為另外，根據(jù)實(shí)際發(fā)酵過(guò)程，系統(tǒng)（8）中向量v和x分別滿足條件：

下面給出系統(tǒng)（8）及其解的一些性質(zhì).

性質(zhì)1 對(duì)任意的k∈K和v∈V，函數(shù)f（x（t），v，k）滿足f∈C（［0，T］，R8）且f在上關(guān)于x是局部Lipschitz連續(xù)的.

性質(zhì)2 對(duì)任意的k∈K和v∈V，函數(shù)f（x（t），v，k）滿足線性增長(zhǎng)條件，即存在常數(shù)α，β＞0使得

定理1 對(duì)任意的k∈K和v∈V，系統(tǒng)（8）存在唯一的解，記為x（·；v，k）.另外，x（·；v，k）在K上關(guān)于k是連續(xù)的.

證明 由性質(zhì)1、2和常微分方程理論［12］知此結(jié)論成立. □

給定x0∈W，記系統(tǒng)（8）的解集為

由以上定理和定義，可得如下定理：

定理2 式（13）定義的可行參數(shù)向量集F是緊集.

2 參數(shù)辨識(shí)模型及優(yōu)化方法

確定系統(tǒng)（8）中的胞內(nèi)動(dòng)力學(xué)參數(shù)對(duì)于明確甘油的代謝機(jī)理、甘油和1，3－PD跨膜運(yùn)輸方式具有非常重要的意義.參數(shù)辨識(shí)模型通常沒(méi)有解析解，因此，構(gòu)造恰當(dāng)?shù)臄?shù)值優(yōu)化方法來(lái)確定參數(shù)辨識(shí)模型的最優(yōu)解是非常必要的.

2.1 參數(shù)辨識(shí)模型

對(duì)甘油歧化連續(xù)發(fā)酵過(guò)程，本文作如下假設(shè).

H1：給定v∈V及任意的x0∈W，k∈F，系統(tǒng)（8）存在唯一的穩(wěn)態(tài).

在假設(shè)H1下，給定v j，j∈Il＝｛1，2，…，l｝，其中l(wèi)表示實(shí)驗(yàn)次數(shù)，已測(cè)得連續(xù)發(fā)酵達(dá)到穩(wěn)態(tài)時(shí)生物量、胞外甘油、胞外1，3－PD、乙醇和乙酸的濃度分別為y j1、y j2、y j3、y j4、y j5.令y j＝（y j1y j2y j3y j4y j5）T∈，j∈Il.另外，當(dāng)系統(tǒng)達(dá)到穩(wěn)態(tài)時(shí)，系統(tǒng)（8）的解x（T j；v j，k）滿足

其中T j表示第j次實(shí)驗(yàn)達(dá)到穩(wěn)態(tài)的時(shí)刻.

在連續(xù)發(fā)酵過(guò)程中，由于堿的加入會(huì)對(duì)乙醇、乙酸的濃度產(chǎn)生影響，本文以生物量、胞外甘油和胞外1，3－PD三種物質(zhì)的計(jì)算值與實(shí)驗(yàn)值之間的平均相對(duì)誤差作為優(yōu)化目標(biāo)，即

這樣，以J（k）作為目標(biāo)函數(shù)，本文建立如下估計(jì)連續(xù)發(fā)酵胞內(nèi)動(dòng)力學(xué)參數(shù)的辨識(shí)模型：

（PIM） minJ（k）

定理3 參數(shù)辨識(shí)模型（PIM）存在最優(yōu)解k＊，即罷k＊∈F有下式成立：

證明 由定理1和2可證此結(jié)論成立. □

2.2 改進(jìn)粒子群算法

粒子群算法（PSO）是由Kennedy和Eberhart于1995年提出的［13］.目前，PSO作為一種新的全局搜索算法已被廣泛用于智能算法、優(yōu)化和其他領(lǐng)域［14～16］.與其他進(jìn)化計(jì)算方法相比，PSO具有搜索能力強(qiáng)、收斂速度快、設(shè)置參數(shù)少、程序?qū)崿F(xiàn)簡(jiǎn)單等特點(diǎn)，既適合科學(xué)研究，又適合工程應(yīng)用.

PSO源于對(duì)鳥(niǎo)群捕食行為的研究.當(dāng)一群鳥(niǎo)在隨機(jī)搜索食物時(shí)，所有的鳥(niǎo)都不知道食物在哪里.但它們知道當(dāng)前位置離食物還有多遠(yuǎn).那么找到食物最簡(jiǎn)單有效的策略就是搜尋目前離食物最近的鳥(niǎo)的周?chē)鷧^(qū)域.PSO就是從這種模型中得到啟發(fā)而產(chǎn)生，并用于解決優(yōu)化問(wèn)題的.PSO最初是處理連續(xù)優(yōu)化問(wèn)題的，目前，其應(yīng)用已擴(kuò)展到組合優(yōu)化問(wèn)題中.在PSO中，每個(gè)優(yōu)化問(wèn)題的解都對(duì)應(yīng)于搜索空間中一只鳥(niǎo)的位置，這些鳥(niǎo)被稱(chēng)之為“粒子”.所有的粒子都有一個(gè)由被優(yōu)化的參數(shù)決定的適應(yīng)值（目標(biāo)函數(shù)），每個(gè)粒子還有一個(gè)速度，決定它們飛翔的方向和距離.在每一次迭代過(guò)程中，粒子通過(guò)跟蹤兩個(gè)最優(yōu)解來(lái)更新自己的位置.一個(gè)是粒子本身所找到的最優(yōu)解，即個(gè)體極值pbest；另一個(gè)是整個(gè)種群目前找到的最優(yōu)解，即全局極值gbest.通過(guò)這種方式，粒子們就追隨當(dāng)前的最優(yōu)粒子在解空間中搜索并朝著好的區(qū)域移動(dòng).

下面給出求解（PIM）的改進(jìn)粒子群算法具體步驟：

式中：分別表示第i個(gè)粒子迭代到第s步時(shí)的速度和位置的第n個(gè)分量，c1、c2是正常數(shù)，r1n、r2n、r3n是［0，1］上服從均勻分布的隨機(jī)數(shù).如果出界，則重新執(zhí)行Step 5操作直到其在界內(nèi)為止.為了避免算法在迭代后期粒子出現(xiàn)震蕩的現(xiàn)象，ω（s）按下式變化：

3 數(shù)值結(jié)果

依實(shí)際發(fā)酵過(guò)程，初始濃度為x01＝0.115 g／L，x02＝495 mmol／L，x03＝x04＝x05＝x06＝x07＝x08＝0.選取連續(xù)發(fā)酵實(shí)驗(yàn)測(cè)得的22組數(shù)據(jù)（其中底物過(guò)量情形的10組和底物限制情形的12組），用改進(jìn)粒子群算法求解（PIM），得到最優(yōu)的胞內(nèi)動(dòng)力學(xué)參數(shù)k＊如表1所示.特別地，本文得到最優(yōu)參數(shù)對(duì)應(yīng)的平均相對(duì)誤差為30.195%，而由參數(shù)k0計(jì)算所得平均相對(duì)誤差為50.830%，從而優(yōu)化過(guò)程使平均相對(duì)誤差減小了20.635%.這里，在改進(jìn)粒子群算法中參數(shù)m、c1、c2、Tmax、ε、ρHJ、ωmax和ωmin分別取值為50、2.0、2.0、500、0.001、0.3、0.9和0.4.另外，將改進(jìn)的粒子群算法與文獻(xiàn)［9］和［11］中算法的收斂性進(jìn)行了比較（如圖1所示）.從圖中可以看出，文獻(xiàn)［9］和［11］中的算法很容易收斂到局部最優(yōu)解，而本文提出的算法具有較好的全局收斂性.

表1 系統(tǒng)（8）中最優(yōu)胞內(nèi)動(dòng)力學(xué)參數(shù)Tab.1 The optimal intracellular kinetic parameters in the System（8）

圖1 3種算法的收斂性曲線Fig.1 Convergence curves of three algorithms

4 結(jié) 論

本文給出了描述連續(xù)發(fā)酵甘油生產(chǎn)1，3－PD過(guò)程的動(dòng)力學(xué)系統(tǒng)，并討論了該系統(tǒng)的一些性質(zhì).為了精確估計(jì)模型中胞內(nèi)動(dòng)力學(xué)參數(shù)，建立了參數(shù)辨識(shí)模型，并證明了最優(yōu)參數(shù)的存在性.最后，構(gòu)造了改進(jìn)粒子群算法求解參數(shù)辨識(shí)模型.數(shù)值結(jié)果表明：計(jì)算值與實(shí)驗(yàn)數(shù)據(jù)的平均相對(duì)誤差減小了20.635%，且該算法比其他算法具有更好的全局收斂性.

［1］修志龍.1，3－丙二醇的微生物法生產(chǎn)分析［J］.現(xiàn)代化工，1999（3）：33－35

［2］BIEBL H，MENZEL K，ZENG An－ping，etal.Microbial production of 1，3－propanediol［J］.Applied Microbiology and Biotechnology，1999，52（3）：289－297

［3］ZENG An－ping，ROSE A，BIEBL H，etal.Multiple product inhibition and growth modeling ofClostridiumbutyricumandKiebsiellapneumoniaein glycerol fermentation［J］.Biotechnology and Bioengineering，1994，44（8）：902－911

［4］修志龍，曾安平，安利佳.甘油生物歧化為1，3－丙二醇過(guò)程的動(dòng)力學(xué)數(shù)學(xué)模擬和多穩(wěn)態(tài)研究［J］.大連理工大學(xué)學(xué)報(bào)，2000，40（4）：428－433

（XIU Zhi－long，ZENG An－ping，AN Li－jia.Mathematical modeling of kinetics and research on multiplicity of glycerol bioconversion to 1，3－propanediol［J］.Journal of Dalian University of Technology，2000，40（4）：428－433）

［5］高彩霞，馮恩民，王宗濤，等.微生物間歇發(fā)酵生產(chǎn)1，3－丙二醇過(guò)程辨識(shí)與優(yōu)化［J］.大連理工大學(xué)學(xué)報(bào)，2006，46（5）：771－774

（GAO Cai－xia，F(xiàn)ENG En－min，WANG Zong－tao，etal.Parameter identification and optimization of process for bio－dissimilation of glycerol to 1，3－propanediol in batch culture［J］.Journal of Dalian University of Technology，2006，46（5）：771－774）

［6］SUN Y Q，QI W T，TENG H，etal.Mathematical modeling of glycerol fermentation byKlebsiella pneumoniae：Concerning enzyme－catalytic reductive pathway and transport of glycerol and 1，3－propanediol across cell membrane［J］.Biochemical Engineering Journal，2008，38（1）：22－32

［7］CHEN Q H，HE G Q，SCHWARZ P.Studies on cultivation kinetics for elastase production byBacillus sp.EL31410［J］.Journal of Agricultural and Food Chemistry，2004，52（11）：3356－3359

［8］FENG Y Y，HE Z M，SONG L F，etal.Kinetics of β－mannanase fermentation by Bacillus licheniformis［J］.Biotechnology Letters，2003，25（14）：1143－1146

［9］TAXARWS A P M，COELHO M A Z，COUTINHO J A P，etal.Laccase improvement in submerged cultivation：induced production and kinetic modeling［J］.Journal of Chemical Technology and Biotechnology，2005，80：669－676

［10］PINCHUK R J，BROWN W A，HUGHES S M，etal.Modeling of biological processes using selfcycling fermentation and genetic algorithms［J］.Biotechnology and Bioengineering，2000，67（1）：19－24

［11］WANG F S，SU T L，JANG H J.Hybrid differential evolution for problems of kinetic parameter estimation and dynamic optimization of an ethanol fermentation process［J］.Industrial ＆Engineering Chemistry Research，2001，40（13）：2876－2885

［12］龐特里亞金ЛC.常微分方程［M］.6版.林武中，倪明康，譯.北京：高等教育出版社，2006

［13］KENNEDY J，EBERHART R.Particle swarm optimization［C］／／Proceedings of the IEEE International Conference on Neural Networks.Piscataway：IEEE，1995：1942－1948

［14］CLERC M，KENNEDY J.The particle swarm：explosion，stability，and convergence in a multidimensional complex space［J］.IEEE Transactions on Evolutionary Computation，2002，6（1）：58－73

［15］YU J B，XI L F，WANG S J.An improved particle swarm optimization for evolving feedforward artificial neural networks［J］.Neural Process Letter，2007，26（3）：217－231

［16］LI H Q，LI L，KIM T H，etal.An improved PSO－based of harmony search for complicated optimization problems［J］.International Journal of Hybrid Information Technology，2008，1（1）：57－64