甘油連續(xù)生物歧化過(guò)程的參數(shù)辨識(shí)與非線性分析

張菁朔,徐恭賢

(渤海大學(xué) 數(shù)學(xué)科學(xué)學(xué)院,遼寧 錦州 121013)

收稿日期：2022-09-06

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目(11101051);遼寧省自然科學(xué)基金資助項(xiàng)目(20180550839);遼寧省教育廳科學(xué)研究項(xiàng)目(LJ2020015)

第一作者：張菁朔(1998—),女,遼寧阜新人,渤海大學(xué)在讀碩士研究生,主要從事最優(yōu)化方法與應(yīng)用研究.

E-mail:3480740062@qq.com

通信作者：徐恭賢(1976—),男,遼寧莊河人,渤海大學(xué)教授,博士,主要從事最優(yōu)化方法與應(yīng)用研究.

E-mail:gxxu@bhu.edu.cn

現(xiàn)如今,甘油生物歧化為1,3-丙二醇過(guò)程的參數(shù)辨識(shí)、最優(yōu)控制等[1-3]生化問(wèn)題受到了學(xué)者們的廣泛關(guān)注.為提高1,3-丙二醇的產(chǎn)率,修志龍等[4]研究了以甘油為原料的連續(xù)發(fā)酵和間歇發(fā)酵過(guò)程,應(yīng)用動(dòng)力學(xué)模型證明了多穩(wěn)態(tài)現(xiàn)象的發(fā)生.王宗濤[5]通過(guò)實(shí)驗(yàn)得到的數(shù)據(jù)及生化反應(yīng)原理,建立了連續(xù)微生物發(fā)酵動(dòng)力學(xué)系統(tǒng)的參數(shù)辨識(shí)模型.李曉紅等[6]研究了微生物間歇發(fā)酵動(dòng)力學(xué)系統(tǒng),采用不可微優(yōu)化理論得到了問(wèn)題最優(yōu)解.宮召華等[7]應(yīng)用離散化方法將所求問(wèn)題轉(zhuǎn)化,并應(yīng)用改進(jìn)的粒子群算法求解參數(shù).王冬妮[8]針對(duì)連續(xù)發(fā)酵生產(chǎn)1,3-丙二醇,建立了分?jǐn)?shù)階微生物發(fā)酵模型以及分?jǐn)?shù)階參數(shù)辨識(shí)優(yōu)化問(wèn)題.文獻(xiàn)[9-11]針對(duì)微生物間歇發(fā)酵過(guò)程的動(dòng)力學(xué)系統(tǒng)建立了不同的參數(shù)辨識(shí)優(yōu)化模型,采用不同的方法求解參數(shù)值.文獻(xiàn)[12]研究了甘油代謝過(guò)程的非線性系統(tǒng)的多目標(biāo)優(yōu)化問(wèn)題.文獻(xiàn)[13]研究了甘油生物合成1,3-丙二醇和2,3-丁二醇的最新工藝,從而使產(chǎn)率最大化.

在對(duì)甘油連續(xù)生物歧化過(guò)程的非線性系統(tǒng)進(jìn)行研究時(shí),本文采用新方法,將甘油連續(xù)生物歧化為1,3-丙二醇過(guò)程的GMA-系統(tǒng)穩(wěn)態(tài)方程作為穩(wěn)態(tài)約束,建立參數(shù)辨識(shí)優(yōu)化模型,應(yīng)用MATLAB序列二次規(guī)劃算法[14]求解模型的最優(yōu)參數(shù)值.將得到的結(jié)果與已有文獻(xiàn)結(jié)果和應(yīng)用LINGO最小二乘法求解的結(jié)果進(jìn)行了比較證明方法的可行性.將辨識(shí)出的參數(shù)值回代到GMA-系統(tǒng)中,計(jì)算了平衡點(diǎn).采用了符號(hào)化實(shí)根隔離法[15]求解出系統(tǒng)所有的平衡點(diǎn).另外,利用線性化方法[16]對(duì)平衡點(diǎn)的穩(wěn)定性進(jìn)行了研究,并討論了不同操作條件下的平衡點(diǎn)的穩(wěn)定性.

1 甘油生物歧化動(dòng)力系統(tǒng)

基于文獻(xiàn)[17],甘油連續(xù)生物歧化為1,3-丙二醇過(guò)程的動(dòng)力學(xué)系統(tǒng)描述如下:

(1)

(2)

(3)

(4)

(5)

其中,t為發(fā)酵時(shí)間,X∈R5,X1為生物量,μ和D分別為細(xì)胞比生長(zhǎng)速率和稀釋速率,CS0和X2分別為初始甘油濃度和反應(yīng)過(guò)程中殘余甘油濃度,X3、X4、X5分別為細(xì)胞外產(chǎn)物1,3-丙二醇、乙酸和乙醇的濃度,qS為底物甘油的比消耗速率,qPD、qHAc、qEtOH分別表示1,3-丙二醇、乙酸和乙醇的比生成速率,μ、qS、qPD、qHAc、qEtOH表達(dá)式如下:

(6)

qS=p1+μ/p2+p3X2/(X2+p4),

(7)

qPD=p5+μp6+p7X2/(X2+p8),

(8)

qHAc=p9+μp10+p11X2/(X2+p12),

(9)

qEtOH=qS[p13/(p14+DX2)+p15/(p16+DX2)],

(10)

令X6=0.28+X2,X7=1-X2/2 039,X8=1-X3/939.5,X9=1-X4/1 026,X10=1-X5/360.9,X11=X2+p4,X12=X2+p8,X13=X2+p12,X14=DX2+p14,X15=DX2+p16,可以得到如下GMA-系統(tǒng)形式:

(11)

(12)

(13)

(14)

(15)

下面定義:

f6(X,D,CS0,p)=0.28+X2,f7(X,D,CS0,p)=1-X2/2 039,

f8(X,D,CS0,p)=1-X3/939.5,f9(X,D,CS0,p)=1-X4/1 026,

f10(X,D,CS0,p)=1-X5/360.9,f11(X,D,CS0,p)=X2+p4,

f12(X,D,CS0,p)=X2+p8,f13(X,D,CS0,p)=X2+p12,

f14(X,D,CS0,p)=DX2+p14,f15(X,D,CS0,p)=DX2+p16.

根據(jù)文獻(xiàn)[17]的實(shí)驗(yàn)結(jié)果,可以得出1組經(jīng)驗(yàn)參數(shù)值,本文將其視為參數(shù)的初始值,具體如下:

p0=(2.2,0.008 2,28.58,11.43,-2.69,67.69,26.59,15.5,-0.97,33.07,

5.74,85.71,0.025,0.06,5.18,50.45)T.

(16)

其中,Pad為R16中的有界閉集,由此可以得到系統(tǒng)(16)的如下性質(zhì):

性質(zhì)1設(shè)p∈Pad,Pad為R16中的有界閉集,則函數(shù)fi(X(t),D,CS0,p)在區(qū)間[t0,tf]上連續(xù)可微,即fi∈C1([t0,tf];R15),且fi關(guān)于p∈Pad連續(xù).

性質(zhì)2若p∈Pad,則非線性動(dòng)力學(xué)系統(tǒng)(16)存在唯一解X(·;p),且X(·;p)在Pad上關(guān)于p連續(xù).

對(duì)于X0,定義系統(tǒng)(16)的解集S(X0)為

S(X0):={X(·;p)∈C([t0,tf];R15)|X(·;p)為系統(tǒng)(16)對(duì)應(yīng)p∈Pad的解}.

(17)

根據(jù)允許集Pad定義,Pad為R16中的有界閉集,從而為緊集.由性質(zhì)1和性質(zhì)2可知,從p∈Pad到X(·;p)∈S(X0)的映射是連續(xù)的.可以得到如下性質(zhì):

性質(zhì)3集合S(X0)是C1([t0,tf];R15)中的緊集.

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

qS(Xj,p)≥0,j=1,2,…,n,qPD(Xj,p)≥0,j=1,2,…,n,qHAc(Xj,p)≥0,j=1,2,…,n,

性質(zhì)4甘油連續(xù)發(fā)酵過(guò)程GMA-系統(tǒng)的參數(shù)辨識(shí)優(yōu)化模型存在最優(yōu)解p*,即存在p*∈Pad,使得對(duì)于任意p∈Pad,都有J(p*)≤J(p).

3 參數(shù)辨識(shí)結(jié)果

本文選取文獻(xiàn)[17]中21組實(shí)驗(yàn)穩(wěn)態(tài)數(shù)據(jù).為求解參數(shù)辨識(shí)優(yōu)化模型中的參數(shù),選用最優(yōu)化算法中的序列二次規(guī)劃算法[14],令Pad={p∈R16|pa∈[0,100],pd∈[-100,0],a=1,2,3,4,6,7,8,10,…,16,d=5,9}.參數(shù)辨識(shí)結(jié)果如表1所示,計(jì)算值誤差對(duì)應(yīng)表2的本文方法誤差.此外,本文還建立了以動(dòng)力學(xué)系統(tǒng)(1)～(5)等號(hào)右端的平方和為目標(biāo)函數(shù)的參數(shù)辨識(shí)優(yōu)化模型進(jìn)行對(duì)照,并應(yīng)用LINGO軟件最小二乘法求解得到的參數(shù)辨識(shí)結(jié)果如表1所示,計(jì)算值誤差對(duì)應(yīng)表2的LINGO計(jì)算誤差.表2為本文方法誤差和已有文獻(xiàn)的誤差的比對(duì)結(jié)果,其中誤差函數(shù)的數(shù)學(xué)定義式如下:

表1 參數(shù)辨識(shí)結(jié)果

表2 誤差比較

(18)

通過(guò)表2的誤差比較可知,應(yīng)用甘油連續(xù)生物歧化過(guò)程的GMA-系統(tǒng),建立以模型計(jì)算值與實(shí)驗(yàn)測(cè)量值誤差最小為目標(biāo)函數(shù)的參數(shù)辨識(shí)優(yōu)化模型產(chǎn)生的濃度誤差小于已有文獻(xiàn)[5,17]的計(jì)算誤差.

4 系統(tǒng)平衡點(diǎn)的數(shù)值計(jì)算

(19)

性質(zhì)5對(duì)任意的v=(CS0,D)∈Vad,系統(tǒng)至少存在1個(gè)平衡點(diǎn).

圖1為GMA-系統(tǒng)中的生物量X1、甘油濃度X2、1,3-丙二醇濃度X3、乙酸濃度X4和乙醇濃度X5在稀釋速率D為0.1 h-1、0.2 h-1和0.3 h-1時(shí)的計(jì)算結(jié)果.當(dāng)稀釋速率D不變時(shí),隨著初始甘油濃度CS0增大,GMA-系統(tǒng)的平衡點(diǎn)個(gè)數(shù)發(fā)生了變化.以稀釋速率D為0.1 h-1的情況為例,隨著CS0逐漸增大,生物量X1、甘油濃度X2、1,3-丙二醇濃度X3、乙酸濃度X4、乙醇濃度X5先是出現(xiàn)了1個(gè)平衡點(diǎn),接著出現(xiàn)了2個(gè)平衡點(diǎn).

a.生物量X1;b.甘油濃度X2;c.1,3-丙二醇濃度X3;d.乙酸濃度X4;e.乙醇濃度X5圖1 不同稀釋速率下的平衡點(diǎn)計(jì)算結(jié)果Fig.1 Calculation results of equilibrium points at different dilution rates

5 系統(tǒng)平衡點(diǎn)的穩(wěn)定性分析

(20)

(21)

(22)

(23)

表3 不同D和CS0情況下的平衡點(diǎn)數(shù)值解及其穩(wěn)定性

6 結(jié)論

針對(duì)甘油連續(xù)歧化為1,3-丙二醇過(guò)程的非線性GMA-系統(tǒng)進(jìn)行了研究,構(gòu)建了包含甘油連續(xù)發(fā)酵過(guò)程GMA-系統(tǒng)的穩(wěn)態(tài)方程為約束的非線性參數(shù)辨識(shí)優(yōu)化問(wèn)題.與已有的參數(shù)估計(jì)方法和應(yīng)用LINGO最小二乘法求出的參數(shù)相比,本文應(yīng)用序列二次規(guī)劃算法得到的最優(yōu)參數(shù)值誤差更小.針對(duì)GMA-系統(tǒng)進(jìn)行了非線性分析,通過(guò)第1步求解出的參數(shù)值,根據(jù)符號(hào)化實(shí)根隔離法求解出不同操作條件下生物系統(tǒng)的全部平衡點(diǎn),給出了不同稀釋速率下的平衡點(diǎn)分布情況.應(yīng)用幾乎線性系統(tǒng)的穩(wěn)定性理論分析了GMA-系統(tǒng)平衡點(diǎn)的穩(wěn)定性,選取了3種不同稀釋速率下的平衡點(diǎn)的數(shù)值進(jìn)行了穩(wěn)定性分析.