甘油連續(xù)生物歧化過(guò)程的狀態(tài)反饋?zhàn)顑?yōu)控制

徐恭賢,王雨竹,李竹茜

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

0 引言

1,3-丙二醇在化工、食品、化妝品和制藥等行業(yè)中具有廣泛的應(yīng)用[1-3]。目前，國(guó)內(nèi)外研究人員主要以羅伊氏乳桿菌、克雷伯氏桿菌、脫氮假單胞菌、丁酸梭菌等歧化甘油的微生物發(fā)酵法來(lái)生產(chǎn)1,3-丙二醇。Menzel等[4]指出，克雷伯氏桿菌歧化甘油生產(chǎn)1,3-丙二醇過(guò)程中存在著明顯的振蕩、抑制等非線性現(xiàn)象。Zeng等[5]，修志龍等[6]和Sun等[7]先后應(yīng)用過(guò)量動(dòng)力學(xué)模型描述了甘油生物歧化過(guò)程中的過(guò)量代謝特性、多態(tài)現(xiàn)象以及酶催化與基因調(diào)控動(dòng)力學(xué)。通過(guò)應(yīng)用代謝工程原理與方法，文獻(xiàn)[8-11]提高了產(chǎn)物1,3-丙二醇的生產(chǎn)強(qiáng)度。針對(duì)生物系統(tǒng)存在的參數(shù)與模型不確定性問(wèn)題，徐恭賢等[12]和Zhu等[13]分別應(yīng)用H∞與μ分析方法，研究了甘油連續(xù)微生物發(fā)酵系統(tǒng)的魯棒控制，取得了較好的控制效果。Xu等[14]和Paranhos等[15]給出了甘油生物歧化過(guò)程過(guò)量動(dòng)力學(xué)系統(tǒng)、酶催化系統(tǒng)以及基因調(diào)控系統(tǒng)的非線性優(yōu)化模型與求解算法，取得了良好的應(yīng)用效果。為了推斷克雷伯氏桿菌歧化甘油過(guò)程的代謝目標(biāo)，徐恭賢等[16]和Xu等[17]提出了多層規(guī)劃模型，并給出了有效的求解算法。應(yīng)用表明，該類方法可以獲得更好的計(jì)算結(jié)果。Yuan等[18]研究了甘油生物歧化系統(tǒng)的最優(yōu)控制問(wèn)題，但該項(xiàng)研究工作沒(méi)有給出控制變量與系統(tǒng)狀態(tài)(或輸出)變量之間的閉環(huán)關(guān)系，因此最優(yōu)控制屬于開(kāi)環(huán)控制。貝泓涵等[19]研究了微生物連續(xù)發(fā)酵生產(chǎn)1,3-丙二醇過(guò)程的線性反饋?zhàn)顑?yōu)控制，但其給出的方法與計(jì)算結(jié)果存在如下問(wèn)題：① 未能給出高質(zhì)量的最優(yōu)解，主要體現(xiàn)在目的產(chǎn)物1,3-丙二醇的濃度最優(yōu)值較低，見(jiàn)本文結(jié)果分析與比較部分的表3；② 設(shè)計(jì)的反饋控制不在約束范圍內(nèi)，見(jiàn)文獻(xiàn)[19]的圖1，原因是在尋找最優(yōu)反饋控制參數(shù)時(shí)沒(méi)有考慮控制變量的有界性。因此有必要給出新的反饋控制模型與方法?；谝陨戏治?，提出適于甘油連續(xù)生物歧化過(guò)程的狀態(tài)反饋?zhàn)顑?yōu)控制模型及其有效求解算法，取得了較好的應(yīng)用效果。

1 甘油連續(xù)生物歧化過(guò)程的數(shù)學(xué)描述

基于文獻(xiàn)[20]，應(yīng)用如下非線性常微分方程系統(tǒng)描述甘油連續(xù)生物歧化過(guò)程：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

t∈[0,T]

(11)

MX(0)=MX0

(12)

MS(0)=MS0

(13)

MPD(0)=MPD0

(14)

MHAc(0)=MHAc0

(15)

MEtOH(0)=MEtOH0

(16)

式(1)—(16)中：MSF為進(jìn)料甘油的濃度，mmol·L-1；d為稀釋速率，h-1；MX為生物量，g·L-1；MS為反應(yīng)器中甘油的濃度，mmol·L-1；MPD為產(chǎn)物1,3-丙二醇的濃度，mmol·L-1；MHAc為產(chǎn)物乙酸的濃度，mmol·L-1；MEtOH為產(chǎn)物乙醇的濃度，mmol·L-1；MX0、MS0、MPD0、MHAc0和MEtOH0為初始濃度；t為時(shí)間，h；T為終端時(shí)刻，h；μX為比生長(zhǎng)速率，h-1；qS為比消耗速率，mmol·g-1·h-1；qPD、qHAc和qEtOH為產(chǎn)物的比生成速率，mmol·g-1·h-1。

2 最優(yōu)控制模型

針對(duì)甘油連續(xù)生物歧化的非線性常微分方程系統(tǒng)(1)—(16)，以最大化終端時(shí)刻t=T時(shí)目的產(chǎn)物1,3-丙二醇濃度MPD(T)為優(yōu)化性能指標(biāo)，以稀釋速率d和進(jìn)料甘油濃度MSF為控制變量，提出了如式(17)—(35)所示的最優(yōu)控制模型。其中，式(29)—(33)是MX、MS、MPD、MHAc和MEtOH的取值范圍；式(34)—(35)是控制變量d和MSF的取值范圍，式(34)的引入是為了防止反應(yīng)過(guò)程發(fā)生“洗出”現(xiàn)象。

maxMPD(T)

(17)

(18)

(19)

(20)

(21)

(22)

t∈[0,T]

(23)

MX(0)=MX0

(24)

MS(0)=MS0

(25)

MPD(0)=MPD0

(26)

MHAc(0)=MHAc0

(27)

MEtOH(0)=MEtOH0

(28)

0.001≤MX≤10

(29)

0≤MS≤2 039

(30)

0≤MPD≤939.5

(31)

0≤MHAc≤1 026

(32)

0≤MEtOH≤360.9

(33)

0.05≤d≤0.5

(34)

100≤MSF≤2 039

(35)

3 狀態(tài)反饋?zhàn)顑?yōu)控制

在甘油連續(xù)生物歧化過(guò)程的非線性常微分方程系統(tǒng)(1)—(16)中，影響目的產(chǎn)物1,3-丙二醇形成的主要因素有2個(gè)，1個(gè)是控制變量稀釋速率d和進(jìn)料甘油濃度MSF，另1個(gè)是生物量MX和反應(yīng)器中甘油濃度MS。因此，設(shè)計(jì)如下?tīng)顟B(tài)反饋控制：

(36)

式(36)中，K∈R2×2為反饋增益矩陣，可表示為

其中，反饋參數(shù)k11、k12、k21和k22的取值范圍為

將式(36)代入式(17)—(35)中得動(dòng)態(tài)優(yōu)化問(wèn)題(記為DOP1)：

maxMPD(T)

(k11MX+k12MS)MEtOH

式(23)—(33)中:

0.05≤k11MX+k12MS≤0.5

100≤k21MX+k22MS≤2039

4 求解方法

為敘述方便，設(shè)

x=(x1,x2,x3,x4,x5)T=

(MX,MS,MPD,MHAc,MEtOH)T

y=(y1,y2,y3,y4)T=

(k11,k12,k21,k22)T

f(x(T),y,T)=x3(T)

g1(x,y)=μXMX-(k11MX+k12MS)MX

g2(x,y)=(k11MX+k12MS)·

(k21MX+k22MS-MS)-qSMX

g3(x,y)=qPDMX-(k11MX+k12MS)MPD

g4(x,y)=qHAcMX-(k11MX+k12MS)MHAc

(k11MX+k12MS)MEtOH

g(x,y)=(g1(x,y),g2(x,y),g3(x,y),

g4(x,y),g5(x,y))T

p1(x,y)=y1x1+y2x2

p2(x,y)=y3x1+y4x2

x0=(MX0,MS0,MPD0,MHAc0,MEtOH0)T

則動(dòng)態(tài)優(yōu)化問(wèn)題DOP1可表示為如下形式(記為DOP2)：

maxf(x(T),y,T)=x3(T)

x(0)=x0

x≤(10,2 039,939.5,1 026,360.9)T

x≥(0.001,0,0,0,0)T

0.05≤p1(x,y)≤0.5

100≤p2(x,y)≤2 039

yl≤y≤yu

t∈[0,T]

由于問(wèn)題DOP2中含有復(fù)雜的狀態(tài)方程dx/dt=g(x,y)，故基于有限元配置形式先將其離散化。

將[0,T]劃分為m個(gè)有限元

[ηi,ηi+1],i=1,2,…,m

其中，0=η1<η2<…<ηm<ηm+1=T。則dx/dt=g(x,y)在有限元[ηi,ηi+1]上的解可近似表示為

其中：xLx(t)為拉格朗日多項(xiàng)式，xLx(tij)=xij；Lx表示[ηi,ηi+1]上狀態(tài)變量x(t)的配置點(diǎn)個(gè)數(shù)。將xLx(t)代入dx/dt=g(x,y)中，可得

(37)

令tik=ηi+Δηiδk，Δηi=ηi-ηi+1，δk∈[0,1]，則式(37)可進(jìn)一步化為

根據(jù)以上推導(dǎo)與分析，將動(dòng)態(tài)優(yōu)化問(wèn)題DOP2轉(zhuǎn)化為如下非線性規(guī)劃問(wèn)題(記為NLP1)：

maxf(x(T),y,T)=x3(T)

i=1,2,…,m；k=1,2,…,Lx

x10=x0

xLx(tij)≤(10,2 039,939.5,1 026,360.9)T

i=1,2,…,m；j=0,1,…,Lx

xLx(tij)≥(0.001,0,0,0,0)T

i=1,2,…,m；j=0,1,…,Lx

0.05≤p1(tij,xij,y)≤0.5

100≤p2(tij,xij,y)≤2 039

yl≤y≤yu

非線性規(guī)劃問(wèn)題NLP1是一個(gè)具有多變量和約束條件的大規(guī)模優(yōu)化問(wèn)題，為了有效求最優(yōu)解，提出如下算法對(duì)其進(jìn)行求解。

算法1：

步驟1給定計(jì)算次數(shù)n。令迭代次數(shù)w=1;

步驟2在第w(w≥1)次迭代，隨機(jī)產(chǎn)生初始可行點(diǎn)，應(yīng)用內(nèi)點(diǎn)法求解非線性規(guī)劃問(wèn)題NLP1，令最優(yōu)解為u(w)，目標(biāo)值為f(w)；

步驟3若w=n，則停止迭代，轉(zhuǎn)步驟4；否則令w=w+1，轉(zhuǎn)步驟2；

采用多初始點(diǎn)策略可以有效解決內(nèi)點(diǎn)法初始點(diǎn)難以確定的問(wèn)題，改善內(nèi)點(diǎn)法的求解效率。

5 結(jié)果與比較

在甘油連續(xù)生物歧化過(guò)程的狀態(tài)反饋控制中，各物質(zhì)的初始濃度分別取為MX0=0.1 g·L-1，MS0=400 mmol·L-1，MPD0=0 mmol·L-1，MHAc0=0 mmol·L-1和MEtOH0=0 mmol·L-1；終端時(shí)刻取為T(mén)=100 h；反饋控制參數(shù)的取值范圍取為

0≤k11≤1 000

0≤k12≤1 000

0≤k21≤1 000

0≤k22≤1 000

通過(guò)算法1可得目的產(chǎn)物1,3-丙二醇濃度的最優(yōu)值為MPD(T)=546.121 499 mmol·L-1，相應(yīng)的最優(yōu)反饋控制參數(shù)如表1所示。表2給出了終端時(shí)刻T=100 h時(shí)各控制變量以及各物質(zhì)濃度的最優(yōu)值。

表1 最優(yōu)反饋控制參數(shù)

表2 T=100 h時(shí)變量的最優(yōu)值

表3為本文算法與文獻(xiàn)[19]方法的計(jì)算結(jié)果比較，其中各物質(zhì)的初始濃度分別取為MX0=0.1 g·L-1，MS0=400 mmol·L-1，MPD0=0 mmol·L-1，MHAc0=0 mmol·L-1和MEtOH0=0 mmol·L-1；終端時(shí)刻取為T(mén)=100 h?？梢?jiàn)，本文獲得的目的產(chǎn)物1,3-丙二醇的濃度為546.121 499 mmol·L-1，是文獻(xiàn)[19]方法結(jié)果的1.279 463倍，說(shuō)明本文方法具有更好的尋優(yōu)性能。

表3 本文方法與文獻(xiàn)[19]方法的計(jì)算結(jié)果

在狀態(tài)反饋控制器的作用下，各控制變量以及各物質(zhì)濃度隨時(shí)間的變化曲線如圖1—7所示，從圖中可以看出，各變量均在給定的約束范圍內(nèi)。

圖1 控制變量(稀釋速率d)曲線

圖2 控制變量(進(jìn)料濃度MSF)曲線

圖3 生物量(MX)曲線

圖4 甘油濃度(MS)曲線

圖5 1,3-丙二醇濃度(MPD)曲線

圖6 乙酸濃度(MHAc)曲線

圖7 乙醇濃度(MEtOH)曲線

6 結(jié)論

研究了甘油連續(xù)生物歧化過(guò)程的狀態(tài)反饋?zhàn)顑?yōu)控制。以最大化終端時(shí)刻目的產(chǎn)物1,3-丙二醇濃度為優(yōu)化性能指標(biāo)，提出了最優(yōu)控制模型。通過(guò)設(shè)計(jì)線性部分狀態(tài)反饋控制器，將最優(yōu)控制模型轉(zhuǎn)化為含反饋控制參數(shù)的動(dòng)態(tài)優(yōu)化問(wèn)題，并提出有效算法對(duì)其進(jìn)行求解，實(shí)現(xiàn)了甘油連續(xù)生物歧化過(guò)程的狀態(tài)反饋?zhàn)顑?yōu)控制。本文算法獲得的產(chǎn)物1,3-丙二醇濃度的最優(yōu)值為546.121 499 mmol·L-1，是已有方法結(jié)果的1.279 463倍，說(shuō)明本文方法具有更好的尋優(yōu)性能。