基于Laplace-Borel 變換的過(guò)程周期性操作的分析

蘇 醒，翟 持

1.衢州學(xué)院，浙江衢州324000；

2.昆明理工大學(xué)化工學(xué)院，云南昆明650500

隨著過(guò)程工業(yè)的深度集成與優(yōu)化，化工過(guò)程的操作與控制呈現(xiàn)高度耦合的趨勢(shì)，反應(yīng)系統(tǒng)的非穩(wěn)態(tài)操作成為最近的研究熱點(diǎn)[1]?；み^(guò)程通常設(shè)計(jì)在定常態(tài)下運(yùn)行，但有研究表明[2]，人為地使操作變量、反應(yīng)物流向和加料位置等因素呈周期性變化，有可能改善反應(yīng)器的時(shí)均性能。從非穩(wěn)態(tài)操作的角度來(lái)看，外界周期操作或者內(nèi)部不穩(wěn)定性有可能變成過(guò)程強(qiáng)化的積極因素，如果對(duì)其加以恰當(dāng)利用，不僅可以減少緩沖罐等單元設(shè)備的投資，還能取得更優(yōu)的工程效果。

輸入變量為參數(shù)的化工過(guò)程擬穩(wěn)態(tài)曲線如圖1所示。該過(guò)程描述發(fā)生在絕熱、全混流反應(yīng)器中的一個(gè)二級(jí)反應(yīng)，模型方程如圖1 右上角所示，其中：q為反應(yīng)物流量，10 m3/h；k為反應(yīng)速率常數(shù)，其值為1.2；V為反應(yīng)器體積，100 m3。反應(yīng)物的進(jìn)料濃度為Cin，kg/m3；縱坐標(biāo)目標(biāo)函數(shù)為出口反應(yīng)物濃度（y），其值越小表示轉(zhuǎn)化率越高。由圖1 可知，周期性地將操作點(diǎn)在A 和B 之間切換使出口平均反應(yīng)物濃度平均值（yavg）小于y，即提高了反應(yīng)物轉(zhuǎn)化率。然而，分析周期性操作時(shí)還需要考慮反應(yīng)的動(dòng)態(tài)變化。Bittanti 等[3]根據(jù)變分原理給出了周期性操作強(qiáng)化過(guò)程的充分條件，并推導(dǎo)得到了π-判據(jù)。Parulekar等[4]運(yùn)用π-判據(jù)研究了多輸入擾動(dòng)情況下的周期操作問(wèn)題。因?yàn)棣?判據(jù)將原系統(tǒng)線性化近似，所以對(duì)高度非線性系統(tǒng)，π-判據(jù)的分析結(jié)果只對(duì)小振幅的周期性輸入有意義[5]。為了獲得高階非線性分析結(jié)果，Kravaris 等[6]采用中心流型理論對(duì)高度非線性過(guò)程的周期性操作進(jìn)行了分析，該方法僅適用于漸進(jìn)穩(wěn)定系統(tǒng)的分析。Zuyev 等[7]使用非線性頻域方法分析周期性操作，實(shí)現(xiàn)了過(guò)程強(qiáng)化的定量分析，但是該分析方法比較復(fù)雜，且對(duì)非周期相應(yīng)分析效果不佳。

圖1 輸入變量為參數(shù)的化工過(guò)程擬穩(wěn)態(tài)曲線Fig.1Steady-state cerve for a process system

本工作以Laplace 變化為基礎(chǔ)，利用Borel 展開分析非線性過(guò)程的周期性操作問(wèn)題，以期有效地分析單擾動(dòng)和多擾動(dòng)周期操作問(wèn)題。基于Laplace-Borel（LB）變換實(shí)現(xiàn)高度非線性系統(tǒng)的高階解析近似[8]，通過(guò)引入積分算子，非線性動(dòng)態(tài)系統(tǒng)可以描述成一個(gè)無(wú)窮級(jí)數(shù)的泛函形式，并且滿足Shuffle操作。鑒于過(guò)程系統(tǒng)引入多個(gè)周期性輸入時(shí)，輸入?yún)?shù)間的關(guān)聯(lián)以及輸入與系統(tǒng)間的耦合使研究系統(tǒng)變得極為復(fù)雜，分別以單輸入和多輸入周期性操作為例，介紹基于LB 變換的周期性操作的分析方法。

1 周期性操作的分析理論

一個(gè)連續(xù)的非線性動(dòng)態(tài)過(guò)程的最優(yōu)周期性問(wèn)題可描述如下：

其中：x∈Rn，是過(guò)程系統(tǒng)的n維狀態(tài)變量；u∈Rm，是系統(tǒng)的m維周期性輸入；φ和ψ分別為系統(tǒng)的等式、不等式約束條件。最優(yōu)周期操作就是尋找到合適的周期性輸入u，使目標(biāo)函數(shù)J最大化。

1.1 非線性傳遞函數(shù)

引入LB 變換以及相應(yīng)的Shuffle 運(yùn)算，將式（1）所述的最優(yōu)周期性操作問(wèn)題映射到復(fù)變函數(shù)域，形成非線性傳遞函數(shù)。當(dāng)系統(tǒng)由多個(gè)微分方程組成，LB 變換可以實(shí)現(xiàn)變量間的分步解耦[9]?，F(xiàn)對(duì)一維的常微分方程進(jìn)行LB 變換，類似地可推廣到多維系統(tǒng)。

其中：f(x)是一個(gè)連續(xù)方程。式（2）可進(jìn)一步寫成通式：

其中：ξ0…ξn為嵌套積分的積分時(shí)間，Ц 為Shuffle 運(yùn)算符，通過(guò)x0-域變換將式（1）轉(zhuǎn)化為代數(shù)運(yùn)算，實(shí)現(xiàn)非線性傳遞函數(shù)構(gòu)建。

與Laplace 變換相似，LB 變換將常微分方程（ODEs）映射到x0域，從而得到代數(shù)方程組，通過(guò)Shuffle運(yùn)算將系統(tǒng)解耦成獨(dú)立的方程，然后逆LB 變換即得到時(shí)域解。由于非線性方程通過(guò)LB 變換得到的是無(wú)窮級(jí)數(shù)，所以運(yùn)用LB 變換得到的解是近似解。常微分方程通過(guò)Taylor 展開可以得到如下形式：

進(jìn)而定義LB 變換，如式（6）所示。常見的初等函數(shù)的LB 表達(dá)式見表1，a和b為常量，F(xiàn)(x0)為象，f(t)為原象。

表1 典型初等函數(shù)LB 變換Table 1LB transform of some functions

1.2 基于LB 變換的泛函展開

當(dāng)系統(tǒng)為多個(gè)ODEs 時(shí)，Picard 迭代積分不能保證收斂，因此，運(yùn)用泛函展開[11]近似求解ODEs時(shí)的收斂問(wèn)題。對(duì)一個(gè)多維系統(tǒng)，如式（7）所示：

其中：x∈Rn，u∈Rm，A∈Rn*n，B∈Rn*m，ξ(x,u)是多項(xiàng)式型的函數(shù)向量。式（7）LB 變換成式（8）：

其中：XL與XN分別為線性和非線性部分的函數(shù)。

由于出現(xiàn)Shuffle 算子，式（8）中的X不能顯式給出，因此將其泛函展開：X=XL+XN(1)+XN(2)...，相應(yīng)地，其非線性部分XN就能逐次展開。與Laplace變換類似，LB 變換能對(duì)非線性系統(tǒng)進(jìn)行頻域分析。本工作重點(diǎn)討論輸入為余弦變化時(shí)，系統(tǒng)的動(dòng)態(tài)響應(yīng)特征，進(jìn)而研究周期性操作問(wèn)題。為方便理解，將泛函展開描述為輸入-輸出型的框圖形式，如圖2所示。其中，XL為式（8）定義的線性傳遞函數(shù)，L 為L(zhǎng)aplace 變換，fi為如下形式：

圖2 輸入與輸出型的泛函展開框圖Fig.2Block diagram of the functional expansion

將式（10）帶入（8）的泛函展開，可以實(shí)現(xiàn)周期性操作的分析。

2 單輸入周期性操作問(wèn)題

2.1 乙醇發(fā)酵-膜分離耦合過(guò)程

工業(yè)上常用運(yùn)動(dòng)發(fā)酵單胞菌（Zymomonasmobilis）或釀酒酵母（Saccharomycescerevisiae）發(fā)酵生產(chǎn)乙醇，但是乙醇是二次代謝產(chǎn)物，隨著發(fā)酵罐內(nèi)產(chǎn)物的累積，微生物生產(chǎn)乙醇過(guò)程會(huì)被抑制，甚至使微生物失活，因此，當(dāng)前的研究熱點(diǎn)是采用原位產(chǎn)物移除技術(shù)及時(shí)移走發(fā)酵過(guò)程產(chǎn)生的乙醇。為研究產(chǎn)物抑制的時(shí)間滯后問(wèn)題，Jobses等[11]建立了結(jié)構(gòu)化模型，將生物質(zhì)分為關(guān)鍵組分K和其余部分G，其中，只有當(dāng)前時(shí)刻系統(tǒng)中K 的含量能決定下一時(shí)刻生物質(zhì)的生成量。結(jié)構(gòu)化模型區(qū)別于傳統(tǒng)的基質(zhì)-生物質(zhì)-產(chǎn)物（s-x-p）模型，需要增加K 關(guān)鍵組分（e），形成如式（11）所示的四維模型。

其中；cs，cx，ce和cp分別為發(fā)酵罐內(nèi)基質(zhì)、生物質(zhì)和關(guān)鍵組分及產(chǎn)物的濃度，kg/m3；D為稀釋速率，h-1；ms為維持微生物生長(zhǎng)常數(shù)，kg/(kg·h)；csf為輸入底物濃度，kg/m3；ks為飽和常數(shù)，其值為0.5 kg/m3；Ysx和Ypx分別是基于基質(zhì)和產(chǎn)品的屈服因子，其值分別為0.024 449 8 kg/kg和0.052 631 5 kg/kg；k1，k2和k3為經(jīng)驗(yàn)常數(shù)，其值分別為16 h-1，0.497 m3/(kg·h)和0.003 83 m6/(kg2·h)。

為消除乙醇抑制作用，根據(jù)滲透膜的選擇性透過(guò)性能，設(shè)計(jì)了如圖3 所示的乙醇發(fā)酵-膜分離耦合發(fā)酵過(guò)程。膜組件右側(cè)使用吹掃蒸汽攜帶走乙醇產(chǎn)品，使發(fā)酵罐內(nèi)的乙醇含量不足以抑制發(fā)酵過(guò)程。發(fā)酵罐轉(zhuǎn)移的乙醇通量為：

圖3 乙醇的膜發(fā)酵過(guò)程Fig.3 The membrane fermentation process

式中：cp和cpM分別為發(fā)酵側(cè)及滲透?jìng)?cè)的產(chǎn)物濃度，kg/m3；F為膜的乙醇滲透性能，為0.128 3 m/h；AM是膜的面積，m2；VF為發(fā)酵側(cè)的體積，為0.003 m3。

圖4 膜生物發(fā)酵過(guò)程。如圖4 所示，由于存在膜的滲透作用，進(jìn)入發(fā)酵側(cè)的稀釋速率D與出發(fā)酵側(cè)的不再相同，分別用Din和Dout表示進(jìn)、出發(fā)酵側(cè)的稀釋速率，h-1；同樣地，用分別用DMin和DMout表示吹掃蒸汽進(jìn)和出滲透?jìng)?cè)的稀釋速率，h-1。乙醇發(fā)酵-膜分離耦合發(fā)酵過(guò)程的動(dòng)態(tài)模型可寫成：

圖4 膜生物發(fā)酵過(guò)程Fig.4The membrane fermentation process

式中：VM為滲透?jìng)?cè)的體積，為0.000 3 m3；ρ為乙醇的密度，其值為789 kg/m3。

2.2 周期變化吹掃蒸汽稀釋速率的膜分離過(guò)程分析

研究發(fā)現(xiàn)，當(dāng)輸入底物濃度（csf）為140 kg/m3且發(fā)酵側(cè)的輸入稀釋率為0.022 h-1時(shí)，發(fā)酵罐內(nèi)累積的乙醇已經(jīng)對(duì)發(fā)酵過(guò)程產(chǎn)生明顯的抑制作用，并且形成自振蕩。若設(shè)計(jì)膜面積為0.24 m2，吹掃蒸汽稀釋速率為0.5 h-1時(shí)，系統(tǒng)不會(huì)形成自振蕩并且趨向穩(wěn)定。相比于沒有膜分離的過(guò)程，耦合膜分離過(guò)程的乙醇產(chǎn)率提高了8.6%。然而，蒸汽滲透膜存在堵塞和污染等問(wèn)題，通過(guò)周期性改變吹掃蒸汽速率，讓膜表面得到更新，減小膜堵塞的影響。周期操作的目標(biāo)是強(qiáng)化生物質(zhì)的乙醇產(chǎn)率（rp），其數(shù)學(xué)表達(dá)式為：

式中，rp為圖4 所示的兩股出料中乙醇質(zhì)量流率在τ周期的時(shí)間平均值，kg/h。當(dāng)膜組件的AM為0.24m2，DMin為0.5 h-1條件下，穩(wěn)態(tài)操作時(shí)的rp為4.590×10-3kg/h。

如果DMin引入余弦變化：δDMin=Amcos(ωt+φ)，由于是單變量周期性操作，可忽略相位角（φ）的影響，將其LB 變換后得到式（10）。結(jié)合式（11）和式（12）建立發(fā)酵-膜分離耦合模型，LB 變換后得到式（15）：

其中：BL[x]∈Rn*n，BL[u]∈Rn，為多項(xiàng)式B[x,u]的線性部分；x(0)=0，是穩(wěn)態(tài)工作點(diǎn)為原點(diǎn)的初值；A0是多項(xiàng)式A[x]的雅克比矩陣。

根據(jù)式（15）可以獲得x0域動(dòng)態(tài)過(guò)程的泛函展開：

其中線性部分（零階泛函展開）可寫成：

式中：τz和τp分別是LB 變換的零點(diǎn)和極點(diǎn)。根據(jù)圖2 進(jìn)行Shuffle 運(yùn)算，可得到高階的泛函展開。

圖5 給出過(guò)程的一階泛函展開，其中自變量為方程式（10）所述的擾動(dòng)頻率和振幅。由圖5 可知，通過(guò)對(duì)吹掃蒸汽DMin進(jìn)行余弦變化，該發(fā)酵-膜分離系統(tǒng)的乙醇產(chǎn)率得到增強(qiáng)，這種增強(qiáng)效果隨著振幅的增加而增加，而隨著頻率的變化，一階偏差影響不大。因此，可以采用工程上容易實(shí)現(xiàn)的變頻方式周期性改變DMin，應(yīng)盡可能增加周期操作的振幅，來(lái)實(shí)現(xiàn)對(duì)該過(guò)程的強(qiáng)化。

圖5 乙醇在發(fā)酵側(cè)的偏移（一階泛函展開）Fig.5The offset of the ethanol in the fermentation side

3 多輸入周期性操作問(wèn)題

3.1 串聯(lián)-平行反應(yīng)過(guò)程

基于LB 變換可以方便地實(shí)現(xiàn)多輸入擾動(dòng)的分析。區(qū)別于單輸入周期性操作，多輸入時(shí)，相位角偏差（φi-φj）是一個(gè)重要的操作變量，因此，多輸入周期性操作問(wèn)題的分析因素較多。在串聯(lián)-平行反應(yīng)中，關(guān)注的核心是中間產(chǎn)物的產(chǎn)量最大化，同時(shí)產(chǎn)生盡量少的副產(chǎn)物?？紤]一個(gè)在絕熱的連續(xù)全混流反應(yīng)器（CSTR）中發(fā)生的反應(yīng)：

其中：k1 和k2 分別為兩個(gè)基元反應(yīng)的反應(yīng)速率常數(shù)。可以建立系統(tǒng)的動(dòng)態(tài)方程：

其中：V為反應(yīng)體積，m3；CAF和CBF為A 和B 的進(jìn)料濃度kmol/ m3；qA和qB為原料A 和B 的體積流率，m3/h。該體系的目標(biāo)是產(chǎn)物P 的時(shí)均產(chǎn)量最大化，同時(shí)考慮副產(chǎn)物R 分離成本及原料A和B 的成本：

式中：J為反應(yīng)體系的總收益的評(píng)價(jià)指標(biāo)；w為副產(chǎn)物R 分離成本相對(duì)于目標(biāo)產(chǎn)物權(quán)重；w1與w2分別為原料A 和B 的成本相對(duì)于目標(biāo)產(chǎn)物的權(quán)重。

3.2 穩(wěn)態(tài)操作

當(dāng)系統(tǒng)穩(wěn)態(tài)時(shí)，目標(biāo)產(chǎn)物最大化僅與兩個(gè)反應(yīng)進(jìn)度（ηi）有關(guān)：

式中：J0為穩(wěn)態(tài)時(shí)產(chǎn)物P 產(chǎn)量最大化函數(shù)；FA0為反應(yīng)物A 的初始摩爾流率，kmol/h；θB為B 與A 摩爾流率之比。通過(guò)解析分析得到各組分的穩(wěn)態(tài)濃度如下：

假設(shè)式中w1，w2，F(xiàn)AF及FBF恒定，令x=η1FA0，可推出J0最大條件及狀態(tài)（?J0/?x=0，?2J0/?x2<0）。

當(dāng)不考慮副產(chǎn)物R 的分離成本，w為0，并且兩個(gè)反應(yīng)常數(shù)之比（ζ=k2/k1）為0.1 時(shí)，可得xmax為0.759 747，ym=η2FA0，其值為0.182 532。取1/θB為0.82、稀釋率（dA=qA/V）為0.1 h-1，k1為1。對(duì)應(yīng)的各物質(zhì)的穩(wěn)態(tài)解為：CA0為0.327 0 kmol/m3，CB0為0.123 8 kmol/m3，CP0為0.208 3 kmol/m3。

3.3 多輸入的周期性操作

現(xiàn)討論反應(yīng)物A和B的濃度周期性變化對(duì)系統(tǒng)的影響。余弦擾動(dòng)U1和U2的LB 變換分別為：

其中：ω1和ω2為兩個(gè)擾動(dòng)的角速度，rad/h；θ為兩個(gè)擾動(dòng)U1和U2的相位移差，h-1；Am1和Am2為振幅，kmol/m3。為方便分析，?。害?=ω2=ω，將式（19）進(jìn)行LB 變換，并帶入（24）得：

通過(guò)變量解耦及消去瞬態(tài)項(xiàng)得到線性項(xiàng)的解滿足如下形式：

由以上參數(shù)得到：τp1=0.7150，τp2=0.3139，τp3=0.1976。τXi分別為不同變量的一階項(xiàng)零點(diǎn)。由此可得一階余項(xiàng)為：A

由上式可知，系統(tǒng)泛函序列取一階非線性展開會(huì)生成余項(xiàng)，發(fā)生對(duì)稱破裂分岔。該系統(tǒng)非線性部分是無(wú)記憶算子，運(yùn)用LB 轉(zhuǎn)換的泛函序列具有收斂性質(zhì)，截取的泛函展開項(xiàng)越多，計(jì)算結(jié)果就越趨近真實(shí)解。

當(dāng)不考慮副產(chǎn)物分離，即w為0，周期性操作僅改變式（19）產(chǎn)物P 的時(shí)均值，兩個(gè)角位移對(duì)P的時(shí)均值的影響如圖6 所示。

圖6 產(chǎn)物P 的一階余項(xiàng)隨頻率和振幅變化關(guān)系Fig.6The first order offset of product Pconcentration with varying forcing amplitude and frequency

由圖6 可知，當(dāng)系統(tǒng)處于最優(yōu)穩(wěn)態(tài)運(yùn)行附近，對(duì)反應(yīng)物A 和B 的濃度實(shí)施周期性操作會(huì)導(dǎo)致產(chǎn)物P 的濃度減小，性能函數(shù)隨著振幅的增加、頻率的降低而惡化，并且θ為0 時(shí)的惡化程度高過(guò)θ為π/2 時(shí)的。該案例也說(shuō)明，不是所有的非線性系統(tǒng)都適合采用周期性操作的方法實(shí)現(xiàn)過(guò)程強(qiáng)化。

4 結(jié) 論

本工作使用LB 傳遞函數(shù)研究連續(xù)乙醇發(fā)酵-膜分離系統(tǒng)，串聯(lián)-平行反應(yīng)過(guò)程的解析近似解，在此基礎(chǔ)上討論了多波動(dòng)參數(shù)的周期型操作對(duì)系統(tǒng)性能函數(shù)的影響。LB 變換是一個(gè)非線性傳變換，將原問(wèn)題轉(zhuǎn)換成一系列泛函展開，通過(guò)Shuffle 運(yùn)算，實(shí)現(xiàn)系統(tǒng)的分部解耦。本工作提出的方法也可以用來(lái)研究參數(shù)對(duì)非線性過(guò)程敏感性的問(wèn)題，對(duì)后續(xù)過(guò)程控制及操作提供相關(guān)的信息。