貝 泓 涵， 張 立 衛(wèi)， 孫 菁

( 大連理工大學(xué) 數(shù)學(xué)科學(xué)學(xué)院, 遼寧 大連 116024 )

0 引 言

1,3-丙二醇在化妝品、聚合物、膠黏劑、潤(rùn)滑劑、藥物等方面有廣泛的應(yīng)用[1]．目前，生產(chǎn)1,3-丙二醇的方法主要有兩種：化學(xué)合成和微生物發(fā)酵．相比于化學(xué)合成，微生物發(fā)酵更具有吸引力，因?yàn)檫^(guò)程更容易施行，并且從環(huán)保角度上來(lái)講沒(méi)有毒副產(chǎn)物產(chǎn)生．然而，由于微生物發(fā)酵法生產(chǎn)1,3-丙二醇相比于化學(xué)合成產(chǎn)量較低，目前還不能應(yīng)用于工業(yè)生產(chǎn)上．這使得越來(lái)越多的學(xué)者聚焦于提高微生物發(fā)酵生產(chǎn)1,3-丙二醇的產(chǎn)量上[2-4]．

自20世紀(jì)80年代起就開(kāi)展了對(duì)1,3-丙二醇的微生物轉(zhuǎn)化過(guò)程研究，近些年學(xué)者們對(duì)該過(guò)程的研究更為細(xì)致精確[5]．微生物發(fā)酵生產(chǎn)1,3-丙二醇主要有3種方式：間歇發(fā)酵、連續(xù)發(fā)酵和批式流加發(fā)酵．連續(xù)發(fā)酵是指在裝有一定體積發(fā)酵液(含有菌種和培養(yǎng)基)的發(fā)酵罐中，以一定速率注入具有一定濃度底物的同時(shí)以相同速率從發(fā)酵罐中取出發(fā)酵液，整個(gè)過(guò)程中保持發(fā)酵液體積不變．其具有生產(chǎn)強(qiáng)度高、產(chǎn)量穩(wěn)定、自動(dòng)化程度高的優(yōu)勢(shì)，因此，本文將著眼于研究連續(xù)發(fā)酵過(guò)程．

對(duì)連續(xù)發(fā)酵過(guò)程的研究學(xué)者們已取得了一定成績(jī)．Zhang等[6]根據(jù)1,3-丙二醇跨細(xì)胞膜的可能傳輸機(jī)制提出了一種描述細(xì)胞內(nèi)還原途徑的非線性混雜系統(tǒng)．Gao等[7]認(rèn)為，甘油通過(guò)細(xì)胞膜有3種可能的途徑，其中包括1,3-丙二醇的被動(dòng)擴(kuò)散和主動(dòng)轉(zhuǎn)運(yùn)，并建立了一種改進(jìn)的14維非線性混合動(dòng)力系統(tǒng)，用遺傳調(diào)控來(lái)描述微生物的連續(xù)發(fā)酵過(guò)程．Lv等[8]在微生物連續(xù)培養(yǎng)過(guò)程中考慮了一種非線性不可微動(dòng)態(tài)系統(tǒng)，包括3-羥基丙醛抑制機(jī)制對(duì)細(xì)胞生長(zhǎng)所有可能的代謝通路，以及甘油和1,3-丙二醇在細(xì)胞膜上的輸運(yùn)系統(tǒng)．本文將研究微生物連續(xù)發(fā)酵過(guò)程的反饋?zhàn)顑?yōu)控制問(wèn)題．

近年來(lái)，學(xué)者們提出了幾種新的求解連續(xù)不等式約束最優(yōu)控制問(wèn)題的計(jì)算方法，包括非光滑牛頓法[9]、約束轉(zhuǎn)錄法[10]、精確罰方法[11]、無(wú)損卷積法和基于區(qū)間分析的方法．雖然這些方法在解決實(shí)際問(wèn)題上已被證明非常有效，但它們只能夠?qū)崿F(xiàn)開(kāi)環(huán)控制，這種開(kāi)環(huán)控制在實(shí)踐中不能保證魯棒性．因此，在微生物連續(xù)發(fā)酵生產(chǎn)1,3-丙二醇的最優(yōu)化問(wèn)題上建立一種反饋控制策略，可以更好地保證魯棒性，并實(shí)現(xiàn)閉環(huán)控制．

解決最優(yōu)反饋控制的傳統(tǒng)方法包括解決HJB偏微分方程(即使對(duì)于無(wú)約束的問(wèn)題來(lái)說(shuō)也是一項(xiàng)困難的任務(wù))．因此，一般而言，計(jì)算約束系統(tǒng)的最優(yōu)反饋控制是非常具有挑戰(zhàn)性的．盡管如此，文獻(xiàn)中還是有部分嘗試．在文獻(xiàn)[12]中，通過(guò)求解一個(gè)反饋增益為決策變量的優(yōu)化問(wèn)題得到了一個(gè)最優(yōu)PID控制器．這個(gè)問(wèn)題包含了由于控制器的工程規(guī)范而產(chǎn)生的連續(xù)的不等式約束，比如超調(diào)和上升時(shí)間的限制．利用約束轉(zhuǎn)錄方法文獻(xiàn)[13]確定了最優(yōu)反饋增益，這是一種處理狀態(tài)約束的著名計(jì)算方法．解決反饋控制問(wèn)題的其他相關(guān)方法包括：計(jì)算魯棒次優(yōu)控制器的敏感罰方法[14-15]，以及鄰近的極值方法[16-17]——當(dāng)系統(tǒng)的新信息可用時(shí)，調(diào)用規(guī)則來(lái)更新最優(yōu)控制．文獻(xiàn)[18]描述了一種專(zhuān)門(mén)為登月模塊設(shè)計(jì)的反饋控制方法，該方法基于時(shí)間尺度和控制參數(shù)化技術(shù)的新組合[19]、著名的彭特里亞金的最小原理[20]和三次樣條逼近[21]．文獻(xiàn)[22]考慮了一個(gè)一般的最優(yōu)控制問(wèn)題，即給定結(jié)構(gòu)的反饋控制器通過(guò)改變某些可調(diào)參數(shù)來(lái)優(yōu)化．

本文提出一種微生物連續(xù)發(fā)酵生產(chǎn)1,3-丙二醇的線性反饋控制策略，建立線性反饋控制器使得1,3-丙二醇的產(chǎn)量最大化．

1 問(wèn)題描述

根據(jù)發(fā)酵實(shí)驗(yàn)過(guò)程，做出如下假設(shè)[23]：

(H1) 發(fā)酵罐中物質(zhì)組成不隨空間位置變化，即反應(yīng)器內(nèi)各點(diǎn)的物質(zhì)組成均一；

(H2) 連續(xù)加入的培養(yǎng)基只含有甘油，反應(yīng)器中的培養(yǎng)液以稀釋速率D輸出．

在上述假設(shè)下，微生物連續(xù)發(fā)酵過(guò)程中生物質(zhì)、基質(zhì)和產(chǎn)物的質(zhì)量平衡關(guān)系可以表示為如下的非線性動(dòng)態(tài)系統(tǒng)：

(1)

并且

xi(0)=x0i；i=1，2，3，4，5

(2)

式中：x1(t)、x2(t)、x3(t)、x4(t)和x5(t)分別表示在時(shí)刻t生物質(zhì)、細(xì)胞外甘油、細(xì)胞外1,3-丙二醇、乙酸和乙醇的濃度(mmol/L)，x0i是其初始濃度(mmol/L)；t∈[0，tf]，tf是終端時(shí)刻；μ是細(xì)胞的比生長(zhǎng)速率(h-1)；q2是底物的比消耗速率(h-1)；q3、q4、q5分別表示1,3-丙二醇、乙酸和乙醇的比生長(zhǎng)速率(h-1)；Cs0表示注入的甘油濃度．

(3)

(4)

(5)

(6)

(7)

其中μm=0.67 h-1，是最大比生長(zhǎng)速率；ks=0.28 mmol/L，是Monod飽和常數(shù)．

令x(t)=(x1(t)x2(t) …x5(t))T，x0=(x01x02…x05)T，u(t)∶=(DCs0)，并且，f(x(t)，u(t))∶=(f1(t) …f5(t))T．

那么非線性控制系統(tǒng)為

x.

x(0)=x0

(8)

考慮到實(shí)際的發(fā)酵過(guò)程中，各物質(zhì)均有臨界濃度，即狀態(tài)向量是有嚴(yán)格區(qū)間的，因此，給出集合W來(lái)界定生物質(zhì)、甘油以及產(chǎn)物的濃度區(qū)間：

(9)

式(9)可以如下式所示等價(jià)地轉(zhuǎn)換為連續(xù)狀態(tài)不等式約束：

hi(x(t)，u(t))≤0；t∈[0，tf]，i=1，…，10

(10)

稀釋速率D以及注入的甘油濃度Cs0有相應(yīng)的臨界值，因此控制變量u(t)應(yīng)滿足如下約束：

u*≤u(t)≤u*；t∈[0，tf]

其中u*和u*分別為u(t)的下確界、上確界．

記x3(·|u(t))為方程(8)的解(1,3-丙二醇濃度)，那么最優(yōu)化問(wèn)題如下所示：

問(wèn)題P0選擇u(t)使得成本函數(shù)最小化

minJ0(u)=-x3(tf|u(t))s.t.

x(t)∈W

u*≤u(t)≤u*

2 反饋控制

在甘油歧化過(guò)程中，對(duì)1,3-丙二醇最終濃度影響最大的兩個(gè)因素就是生物質(zhì)和甘油的濃度，并且線性狀態(tài)反饋也是常見(jiàn)的反饋控制結(jié)構(gòu)之一．因此，選用與生物質(zhì)、甘油濃度相關(guān)的線性反饋控制器：

u(t)=(DCs0)=(φ1(x(t)，ξ)φ2(x(t)，ξ))

(11)

其中φ1(x(t)，ξ)=ξ1x1(t)+ξ2x2(t)，φ2(x(t)，ξ)=ξ3x1(t)+ξ4x2(t)；記ξ=(ξ1ξ2ξ3ξ4)T∈R4是反饋控制參數(shù)向量；φ1和φ2是如上的連續(xù)可微函數(shù)．那么，反饋參數(shù)ξ1、ξ2、ξ3、ξ4是需要最優(yōu)選擇的決策變量．

反饋參數(shù)須滿足如下有界約束條件：

ξ=(ξ1ξ2ξ3ξ4)T∈Uad=

[α1，β1]×[α2，β2]×[α3，β3]×[α4，β4]

(12)

將式(11)代入式(8)中有

(13)

其中

f～

q2x2(t)

(14)

方程(14)應(yīng)滿足方程(8)的初始條件，另記x(·|ξ)為方程(14)的解．此外，約束條件(10)變?yōu)?/p>

hi(x(t)|ξ)≤0；t∈[0，tf]，i=1，…，10

(15)

本文的目標(biāo)是通過(guò)狀態(tài)反饋策略使1,3-丙二醇的最終濃度最大化，那么首先考慮在滿足約束條件(15)的基礎(chǔ)上，選擇最優(yōu)反饋參數(shù)ξk(k=1，2，3，4)使得系統(tǒng)成本最小化．

問(wèn)題P 選擇ξ∈Uad使得成本函數(shù)最小化

minJ1(ξ)=-x3(tf|ξ)s.t.x(t)=f～(x(t),ξ)x(0)=x0ξ∈Uad

3 精確罰方法

問(wèn)題P是一個(gè)非線性?xún)?yōu)化問(wèn)題，在此問(wèn)題中，有限數(shù)量的決策變量(反饋控制參數(shù))需要滿足約束條件并進(jìn)行優(yōu)化．這是個(gè)非常難以解決的優(yōu)化問(wèn)題，因?yàn)槊總€(gè)連續(xù)不等式約束條件實(shí)際上都是由無(wú)限個(gè)(對(duì)應(yīng)時(shí)間區(qū)間[0，tf]的每個(gè)點(diǎn))約束條件組成的．那么，問(wèn)題P就可以看成是一個(gè)半無(wú)限優(yōu)化問(wèn)題．為解決這個(gè)問(wèn)題，接下來(lái)使用精確罰方法[24]得到問(wèn)題P的近似問(wèn)題．

x(t)∈W等價(jià)于G(ξ)=0，其中

(16)

顯然G(ξ)=0當(dāng)且僅當(dāng)x(t)∈W．然而，等式約束條件(16)當(dāng)hi=0時(shí)非光滑，因此，一般的優(yōu)化程序在處理此類(lèi)等式約束問(wèn)題上有一定的困難．那么，考慮如下的精確罰函數(shù)：

Jσ(ξ，

)=J1(ξ);=0,G(ξ)=0J1(ξ)+-αG(ξ)+σβ;∈(0,-]∞;=0, G(ξ)≠0ì?í????

(17)

罰函數(shù)Jσ用于懲罰G(ξ)中的較大值，因此，最小化Jσ能讓可行點(diǎn)滿足約束條件(15)．在此基礎(chǔ)上，問(wèn)題P可近似于如下罰問(wèn)題：

問(wèn)題Q 選擇(ξ，)∈Uad×[0，]使得罰函數(shù)Jσ(ξ，)最小化．

minJσ(ξ，)=J1(ξ)+-αG(ξ)+σβs.t.x(t)=(x(t)，ξ)

x(0)=x0

ξ(t)∈Uad

解決這個(gè)問(wèn)題的數(shù)值算法利用了成本函數(shù)的梯度來(lái)找到滿足搜索空間可行區(qū)域的上升方向，采用如下方法來(lái)計(jì)算成本函數(shù)的梯度：

首先，對(duì)于每個(gè)k=1，…，4，考慮如下的變分方程：

?.k=?f～(x(t),ξ)?x?k(t)+?f～(x(t),ξ)?ξk

(18)

φk(0)=0

(19)

其中

?f～(x(t),ξ)?x=

(20)

并且

?f～(x(t),ξ)?ξk=?f(x(t),φ(x(t),ξ))?u?φ(x(t),ξ)?ξk

(21)

記φk(·|ξ)為方程(18)、(19)的解，有如下定理：

定理1對(duì)每一ξ∈Γ

(22)

基于定理1，使用標(biāo)準(zhǔn)微分法則，可以得到Jσ的偏微分．

定理2問(wèn)題Q中，Jσ的偏微分如下：

(23)

(24)

其中

(25)

(26)

以及，Φ(x(tf|ξ))=-x3(tf|ξ)．

4 算法及數(shù)值結(jié)果

基于上述內(nèi)容，給出求解問(wèn)題Q的算法：

步驟2令

步驟3以(ξ0，0)作為初值，使用非線性?xún)?yōu)化算法二次序列規(guī)劃方法(SQP)求解問(wèn)題Q．記(ξ*，*)為求得的局部最小解，其中

步驟4如果*<ρ，則停止，取ξ*作為問(wèn)題Q的一個(gè)解．否則，令10σ→σ轉(zhuǎn)到步驟5．

步驟5如果σ≤σmax，那么令(ξ*，*)→(ξ0，0)并轉(zhuǎn)到步驟2．否則停止——t算法無(wú)法得到問(wèn)題Q的解．

在37 ℃，pH=7．0的厭氧條件下，系統(tǒng)的取值如表1所示[25]．

表1 參數(shù)取值

在微生物發(fā)酵過(guò)程控制中，取狀態(tài)向量的臨界值為x*=(0.001 100 0 0 0)T，x*=(10

2 039 939.5 1 026 360.9)；生物質(zhì)、甘油以及1,3-丙二醇、乙酸、乙醇的初始濃度分別為x01=0.1 mmol/L，x02=400 mmol/L，x03=0，x04=0，x05=0．控制變量D、Cs0應(yīng)滿足D(t)∈[0.05，0.67]，Cs0(t)∈[100，1 800]．連續(xù)發(fā)酵過(guò)程中基質(zhì)充分，取反應(yīng)時(shí)間為100 h.在數(shù)值實(shí)驗(yàn)過(guò)程中，作為經(jīng)驗(yàn)數(shù)值，ξ的參數(shù)取值范圍為[0，400]×[0，2]×[0，20]×[0，0.000 1]．通過(guò)本文的算法計(jì)算得出，反饋控制最優(yōu)參數(shù)ξ*=(394.132 8 0.630 4 0.085 0.000 05)，1,3-丙二醇在終端時(shí)刻的濃度為426.836 4 mmol/L．控制變量D、Cs0如圖1、2所示；生物質(zhì)、甘油、1,3-丙二醇、乙酸、乙醇的濃度如圖3所示．

圖1 控制變量DFig.1 Control variable D

圖2 控制變量Cs0Fig.2 Control variable Cs0

(a) 生物質(zhì)

(b) 甘油

(c) 1,3-丙二醇

(d) 乙酸

(e) 乙醇

圖3 生物質(zhì)、甘油、1,3-丙二醇、乙酸、乙醇的濃度變化

Fig.3 The concentration change of biomass, glycerol, 1,3-PD, acetate and ethanol

5 結(jié) 語(yǔ)

本文提出了在微生物發(fā)酵生產(chǎn)1,3-丙二醇過(guò)程中的一種線性反饋控制策略．該策略可以很好地保證魯棒性，并通過(guò)精確罰方法給出了原問(wèn)題的近似問(wèn)題從而使用一種非線性?xún)?yōu)化算法二次序列規(guī)劃方法求得全局最優(yōu)解，而且通過(guò)線性反饋策略?xún)?yōu)先實(shí)現(xiàn)了閉環(huán)控制．希望未來(lái)可以實(shí)現(xiàn)多階段的反饋控制，令反饋參數(shù)為與時(shí)間t相關(guān)的函數(shù)，從而更加精確地實(shí)現(xiàn)控制．

[1] LIU Chongyang, GONG Zhaohua, FENG Enmin,etal. Modelling and optimal control for nonlinear multistage dynamical system of microbial fed-batch culture [J].JournalofIndustrialandManagementOptimization, 2009,5(4):835-850.

[2] LIU Chongyang, GONG Zhaohua, TEO K L. Robust parameter estimation for nonlinear multistage time-delay systems with noisy measurement data [J].AppliedMathematicalModelling, 2018,53:353-368.

[3] YUAN Jinlong, ZHANG Yuduo, YE Jianxiong,etal. Robust parameter identification using parallel global optimization for a batch nonlinear enzyme-catalytic time-delayed process presenting metabolic discontinuities [J].AppliedMathematicalModelling, 2017,46:554-571.

[4] YUAN Jinlong, ZHANG Xu, ZHU Xi,etal. Modelling and pathway identification involving the transport mechanism of a complex metabolic system in batch culture [J].CommunicationsinNonlinearScienceandNumericalSimulation, 2014,19(6):2088-2103.

[5] 修志龍. 微生物發(fā)酵法生產(chǎn)1,3-丙二醇的研究進(jìn)展[J]. 微生物學(xué)通報(bào), 2000,27(4):300-302.

XIU Zhilong. Progress in the production of 1,3-PD by microorganism fermentation [J].MicrobiologyChina, 2000,27(4):300-302. (in Chinese)

[6] ZHANG Yuduo, FENG Enmin, XIU Zhilong. Robust analysis of hybrid dynamical systems for 1,3-propanediol transport mechanisms in microbial continuous fermentation [J].MathematicalandComputerModelling, 2011,54(11):3164-3171.

[7] GAO Kuikui, ZHANG Xu, FENG Enmin,etal. Sensitivity analysis and parameter identification of nonlinear hybrid systems for glycerol transport mechanisms in continuous culture [J].JournalofTheoreticalBiology, 2014,347:137-143.

[8] LV Jiajia, PANG Liping, FENG Enmin. Asymptotical stability of a nonlinear non-differentiable dynamic system in microbial continuous cultures [J].AdvancesinDifferenceEquations, 2017,2017:256. https://doi.org/10.1186/s13662-017-1288-X.

[9] GERDTS M, KUNKEL M. A nonsmooth Newton′s method for discretized optimal control problems with state and control constraints [J].JournalofIndustrialandManagementOptimization, 2008,4(2):247-270.

[10] WANG L Y, GUI W H, TEO K L,etal. Time delayed optimal control problems with multiple characteristic time points: computation and industrial applications [J].JournalofIndustrialandManagementOptimization, 2009,5(4):705-718.

[11] YU Changjun, LI Bin, LOXTON R,etal. Optimal discrete-valued control computation [J].JournalofGlobalOptimization, 2013,56:503-518.

[12] LI Bin, TEO K L, LIM C C,etal. An optimal PID controller design for nonlinear constrained optimal control problems [J].DiscreteandContinuousDynamicalSystems—SeriesB, 2011,16:1101-1117.

[13] LOXTON R, TEO K L, REHBOCK V,etal. Optimal control problems with a continuous inequality constraint on the state and the control [J].Automatica, 2009,45:2250-2257.

[14] LOXTON R, TEO K L, REHBOCK V. Robust suboptimal control of nonlinear systems [J].AppliedMathematicsandComputation, 2011,217:6566-6576.

[15] REHBOCK V, TEO K L, JENNINGS L S. A computational procedure for suboptimal robust controls [J].DynamicsandControl, 1992,2(4):331-348.

[16] FISHER M E, GRANTHAM W J, TEO K L. Neighbouring extremals for nonlinear systems with control constraints [J].DynamicsandControl, 1995,5:225-240.

[17] JIANG C, TEO K L, LOXTON R,etal. A neighboring extremal solution for an optimal switched impulsive control problem [J].JournalofIndustrialandManagementOptimization, 2012,8:591-609.

[18] ZHOU Jingyang, TEO K L, ZHOU Di,etal. Nonlinear optimal feedback control for lunar module soft landing [J].JournalofGlobalOptimization, 2012,52(2):211-227.

[19] LOXTON R, TEO K L, REHBOCK V. Optimal control problems with multiple characteristic time points in the objective and constraints [J].Automatica, 2008,44:2923-2929.

[20] PONTRYAGIN L S.TheMathematicalTheoryofOptimalProcesses[M]. New York: Gordan and Breach, 1986.

[21] BURDEN R L, FAIRES J D.NumericalAnalysis[M]. Boston: Brooks Cole, 2010.

[22] LIN Q, LOXTON R, TEO K L,etal. Optimal feedback control for dynamic systems with state constraints: An exact penalty approach [J].OptimizationLetters, 2014,8(4):1535-1551.

[23] GAO C X, WANG Z T, FENG E M,etal. Nonlinear dynamical systems of bio-dissimilation of glycerol to 1, 3-propanediol and their optimal controls [J].JournalofIndustrialandManagementOptimization, 2005,1(3):377-388.

[24] YU C, TEO K L, ZHANG L,etal. A new exact penalty function method for continuous inequality constrained optimization problems [J].JournalofIndustrialandManagementOptimization, 2010,6:895-910.

[25] 修志龍,曾安平,安利佳. 甘油生物歧化過(guò)程動(dòng)力學(xué)數(shù)學(xué)模擬和多穩(wěn)態(tài)研究[J]. 大連理工大學(xué)學(xué)報(bào), 2004,40(4):428-433.

XIU Zhilong, ZENG Anping, AN Lijia. Mathematical modeling of kinetics and research on multiplicity of glycerol bioconversion to 1,3-propanediol [J].JournalofDalianUniversityofTechnology, 2000,40(4):428-433. (in Chinese)