戴吟 劉飛

基金項(xiàng)目：國(guó)家自然科學(xué)基金項(xiàng)目（批準(zhǔn)號(hào)：61833007）資助的課題。

作者簡(jiǎn)介：戴吟（1998-），碩士研究生，從事工業(yè)過程控制及應(yīng)用的研究。

通訊作者：劉飛（1965-），教授，從事先進(jìn)過程控制方面的教學(xué)與科研工作，fliu@jiangnan.edu.cn。

引用本文：戴吟，劉飛.基于穩(wěn)定逆的生化過程操作曲線求取方法[J].化工自動(dòng)化及儀表，2023，50（6）：000-000.

DOI：10.20030/j.cnki.1000-3932.202306000

摘? 要? 針對(duì)生化過程中操作曲線的求取問題，提出一種非線性穩(wěn)定逆的近似求解方法。該方法根據(jù)過程動(dòng)力學(xué)模型，利用改進(jìn)的穩(wěn)定逆算法來求解對(duì)應(yīng)生產(chǎn)輸出目標(biāo)的操作曲線。并針對(duì)穩(wěn)定逆算法計(jì)算量過大的問題，通過對(duì)穩(wěn)定逆原理中的內(nèi)動(dòng)態(tài)進(jìn)行離散化處理，將原本利用非線性系統(tǒng)內(nèi)動(dòng)態(tài)方程求解內(nèi)動(dòng)態(tài)解析解的過程，利用遞推公式轉(zhuǎn)化成求解內(nèi)動(dòng)態(tài)的離散數(shù)值解，從而加快穩(wěn)定逆算法的計(jì)算速度。最后以化學(xué)反應(yīng)器為例進(jìn)行操作曲線設(shè)計(jì)的仿真，驗(yàn)證了所提方法的可行性。

關(guān)鍵詞? 操作曲線? 動(dòng)力學(xué)模型? 非線性穩(wěn)定逆? 內(nèi)動(dòng)態(tài)方程? 化學(xué)反應(yīng)器

中圖分類號(hào)? TP273? ? ? ? ?文獻(xiàn)標(biāo)志碼? A? ? ? ? ?文章編號(hào)? 1000-3932（2023）06-0000-00

生物及化工生產(chǎn)目標(biāo)的實(shí)現(xiàn)取決于過程溫度、流量及液位等關(guān)鍵變量的操作曲線，一旦工藝確定，控制系統(tǒng)即可按照生產(chǎn)操作曲線進(jìn)行實(shí)時(shí)控制，從而完成最終的生產(chǎn)要求[1]。生化過程操作曲線的選取與產(chǎn)品質(zhì)量、產(chǎn)量及生產(chǎn)時(shí)間等目標(biāo)密切相關(guān)，由此過程工藝的制定是實(shí)現(xiàn)智能生產(chǎn)的核心環(huán)節(jié)，也是整個(gè)生化生產(chǎn)的基礎(chǔ)。一般情況，調(diào)控參數(shù)可以通過工藝人員反復(fù)實(shí)驗(yàn)得到[2～4]，但實(shí)驗(yàn)方法需要耗費(fèi)大量的人力、物力和時(shí)間，尤其難以適應(yīng)現(xiàn)代生產(chǎn)對(duì)多品種、小批量、多目標(biāo)的生產(chǎn)需求，同時(shí)也很難通過實(shí)驗(yàn)設(shè)計(jì)的方法來得到動(dòng)態(tài)操作曲線[5]。根據(jù)生產(chǎn)目標(biāo)，在已知生產(chǎn)過程動(dòng)力學(xué)的基礎(chǔ)上，通過解析的方法求解操作曲線，是值得開拓的研究領(lǐng)域。從控制理論的角度，這一問題可轉(zhuǎn)化為在動(dòng)力學(xué)模型基礎(chǔ)上，根據(jù)模型輸出軌跡（工藝目標(biāo)）求取相應(yīng)的模型輸入軌跡（調(diào)控曲線），即模型逆的存在性及求解算法。

1996年HUNT L R等針對(duì)非最小相位過程首次提出非因果穩(wěn)定逆的方法[6]。由于穩(wěn)定逆方法需要未來無窮時(shí)刻的軌跡，因此改進(jìn)的基于預(yù)覽信息的穩(wěn)定逆方法被提出[7，8]，進(jìn)一步引入輸出誤差和控制器輸出能量為性能指標(biāo)，給出了最優(yōu)預(yù)覽穩(wěn)定逆方法[9]?？紤]到穩(wěn)定逆算法需要很大預(yù)作用時(shí)間區(qū)間，利用最優(yōu)狀態(tài)轉(zhuǎn)移的方法使過程狀態(tài)轉(zhuǎn)移到穩(wěn)定逆初態(tài)，進(jìn)而提出了基于最優(yōu)狀態(tài)轉(zhuǎn)移的穩(wěn)定逆算法[10]。

由于穩(wěn)定逆算法所需計(jì)算量很大，近年來有學(xué)者利用迭代學(xué)習(xí)結(jié)合樣條技術(shù)得到精確跟蹤所需的輸入控制量[11]，還有基于最小二乘、B樣條等基函數(shù)的方法來求解非最小相位過程的穩(wěn)定逆[12]。一般而言，與線性穩(wěn)定逆算法有所區(qū)別，非線性過程穩(wěn)定逆計(jì)算中需要利用Picard迭代方法來得到非線性內(nèi)動(dòng)態(tài)方程的解[8]，對(duì)于復(fù)雜過程模型或者輸出軌跡，求解內(nèi)動(dòng)態(tài)方程時(shí)會(huì)遇到計(jì)算量過大，甚至解析解無法描述的困難。

考慮生化過程非線性動(dòng)力學(xué)情形，筆者基于穩(wěn)定逆方法，提出了一種生產(chǎn)操作曲線的快速求取方法，即根據(jù)過程動(dòng)力學(xué)模型，利用改進(jìn)的穩(wěn)定逆算法來求解對(duì)應(yīng)生產(chǎn)輸出目標(biāo)的操作曲線。首先將穩(wěn)定逆算法中內(nèi)動(dòng)態(tài)進(jìn)行離散化處理，利用遞推公式推導(dǎo)內(nèi)動(dòng)態(tài)的離散數(shù)值表達(dá)式；再以Picard迭代進(jìn)行遞推求取非線性內(nèi)動(dòng)態(tài)方程的數(shù)值解，從而減少計(jì)算量、加快穩(wěn)定逆算法的計(jì)算速度；最后以化學(xué)反應(yīng)器為例，進(jìn)行操作曲線設(shè)計(jì)驗(yàn)證。

1? 操作曲線求取問題

考慮一個(gè)化學(xué)反應(yīng)器過程，包括以下串聯(lián)/并聯(lián)反應(yīng)[13]：

（1）

反應(yīng)器的動(dòng)力學(xué)模型包括物料A、B的摩爾平衡和能量平衡：

（2）

其中，CA和CB是物料A和B的摩爾濃度；T是反應(yīng)器溫度；是進(jìn)料流量；TC是冷卻夾套的溫度；Cp和分別是反應(yīng)混合物的比熱容和密度；是反應(yīng)熱。反應(yīng)速率系數(shù)k1～k3與反應(yīng)器溫度和過程參數(shù)有關(guān)，具體過程參數(shù)見表1，反應(yīng)速率系數(shù)具體計(jì)算公式如下：

（3）

令過程模型輸入分別為進(jìn)料流量u1=F和反應(yīng)器冷卻夾套溫度u2=TC，模型輸出分別是物料B的摩爾濃度y1=CB和反應(yīng)器溫度y2=T。反應(yīng)器穩(wěn)態(tài)平衡點(diǎn)為和，此時(shí)對(duì)應(yīng)于和。

化學(xué)反應(yīng)器的生產(chǎn)控制目標(biāo)是反應(yīng)器產(chǎn)物B的濃度，調(diào)控手段是調(diào)節(jié)進(jìn)料流量和反應(yīng)器冷卻夾套溫度。具體生產(chǎn)過程目標(biāo)是在保持反應(yīng)器溫度恒定的情況下，要求物料B的濃度隨生產(chǎn)時(shí)間按照以下設(shè)定曲線變化：

（4）

其中，tf是調(diào)控生產(chǎn)目標(biāo)的終端時(shí)間。

文中操作曲線設(shè)計(jì)即根據(jù)上述生產(chǎn)需求利用穩(wěn)定逆方法來獲取進(jìn)料流量和反應(yīng)器冷卻夾套溫度的操作曲線。

2? 非線性穩(wěn)定逆近似求解的改進(jìn)算法

2.1? 非線性穩(wěn)定逆

考慮如下多輸入多輸出非線性過程：

（5）

其中，狀態(tài)，控制輸入，輸出，、均為光滑向量場(chǎng)即非線性平滑函數(shù)的向量集合，為非線性平滑函數(shù)。

設(shè)為期望的輸出軌跡，需要找出對(duì)應(yīng)的輸入軌跡和狀態(tài)軌跡，即：

（6）

對(duì)式（5）中的函數(shù)h（x）沿向量場(chǎng)g和f進(jìn)行求導(dǎo)運(yùn)算，即求李導(dǎo)數(shù)和[5]。根據(jù)李導(dǎo)數(shù)定義將輸出對(duì)時(shí)間t求導(dǎo)r次直至出現(xiàn)，其中r為過程的相對(duì)階：

（7）

由相對(duì)階定義知是可逆的，因此可以得到輸出控制律形式：

（8）

外動(dòng)態(tài)定義為過程輸出及其各階導(dǎo)數(shù)：

（9）

將過程狀態(tài)向量的其他n-r個(gè)元素定義為內(nèi)動(dòng)態(tài)：

（10）

將稱為一個(gè)局部微分同胚，即通過非線性坐標(biāo)變換將系統(tǒng)狀態(tài)變量x映射到的新坐標(biāo)系下，并且是可逆的，即存在函數(shù)，于是式（5）可變成坐標(biāo)下的動(dòng)態(tài)系統(tǒng)：

（11）

如果滿足輸出跟蹤軌跡，即，則外動(dòng)態(tài)可以寫成：

（12）

在滿足輸出跟蹤和新坐標(biāo)系的條件下，根據(jù)式（8）得到相應(yīng)的控制輸入：

（13）

由于已知，從而穩(wěn)定逆的計(jì)算就是求解內(nèi)部狀態(tài)參考軌跡，即內(nèi)部狀態(tài)的有界解，可以利用非線性內(nèi)動(dòng)態(tài)方程求解內(nèi)動(dòng)態(tài)。

對(duì)于非線性內(nèi)動(dòng)態(tài)方程，可將其寫為：

（14）

（15）

假設(shè)矩陣Q中穩(wěn)定部分和不穩(wěn)定部分已經(jīng)分開，即：

（16）

其中，Qs所有特征值在左半平面，屬于穩(wěn)定部分；Qu所有特征值在右半平面，屬于不穩(wěn)定部分。

定義如下有界狀態(tài)轉(zhuǎn)移矩陣：

（17）

則內(nèi)動(dòng)態(tài)方程（14）等價(jià)于以下積分方程：

（18）

通過Picard迭代方法可以得到該方程的解[7]：

（19）

當(dāng)?shù)螖?shù)時(shí)，。

文獻(xiàn)[8]對(duì)非線性過程提出了預(yù)覽穩(wěn)定逆方法，只需要知道未來一段時(shí)間內(nèi)的輸出軌跡信息，其中預(yù)覽時(shí)間越長(zhǎng)，跟蹤誤差越小。即積分方程（18）可以近似為：

（20）

設(shè)初始時(shí)刻為t0，利用式（19）求解內(nèi)動(dòng)態(tài)在未來時(shí)間的解析解時(shí)，隨著迭代次數(shù)m的增加，結(jié)果逐漸收斂于實(shí)際值，但的解析表達(dá)式會(huì)趨于復(fù)雜甚至難以描述，造成計(jì)算量越來越大直至無法計(jì)算的應(yīng)用問題；另一方面，生化過程復(fù)雜的動(dòng)力學(xué)模型或輸出軌跡也會(huì)使式（18）的積分方程變得更加復(fù)雜，導(dǎo)致Picard迭代方法收斂困難。

2.2? 非線性穩(wěn)定逆的離散化近似求解算法

2.2.1? 內(nèi)動(dòng)態(tài)離散化處理

假設(shè)內(nèi)動(dòng)態(tài)在未來時(shí)間內(nèi)按時(shí)間均勻分解n段，即化成n+1個(gè)離散的數(shù)值向量，在每一段長(zhǎng)度為時(shí)間內(nèi)，被近似認(rèn)為是一個(gè)定值，而不是原先連續(xù)隨時(shí)間變化的函數(shù)，離散化后的內(nèi)動(dòng)態(tài)向量形式表示如下：

（21）

其中，的形式表示函數(shù)在時(shí)間內(nèi)離散后長(zhǎng)度為n+1的向量。

2.2.2? 遞推內(nèi)動(dòng)態(tài)向量

將式（19）積分方程分為穩(wěn)定內(nèi)動(dòng)態(tài)和不穩(wěn)定內(nèi)動(dòng)態(tài)積分方程，設(shè)，分別以特征值為和（）的內(nèi)動(dòng)態(tài)子系統(tǒng)為例推導(dǎo)不穩(wěn)定內(nèi)動(dòng)態(tài)和穩(wěn)定內(nèi)動(dòng)態(tài)：

（22）

其中，下標(biāo)s、u分別表示穩(wěn)定部分和不穩(wěn)定部分；，M為可逆的內(nèi)動(dòng)態(tài)分解矩陣。

不穩(wěn)定內(nèi)動(dòng)態(tài)向量值的遞推公式如下：

（23）

同理，可以推導(dǎo)出穩(wěn)定內(nèi)動(dòng)態(tài)向量值的遞推公式：

（24）

上述對(duì)于不穩(wěn)定內(nèi)動(dòng)態(tài)和穩(wěn)定內(nèi)動(dòng)態(tài)的遞推方向相反，由于和為第m-1次Picard迭代出的數(shù)值向量，當(dāng)求第m次積分迭代時(shí)每一個(gè)遞推公式中的積分求解結(jié)果均是數(shù)值形式，這樣處理避免了求解復(fù)雜的解析解形式，極大減少了計(jì)算量。

當(dāng)Picard迭代次數(shù)m增加到不穩(wěn)定內(nèi)動(dòng)態(tài)和穩(wěn)定內(nèi)動(dòng)態(tài)收斂時(shí)，最后的輸出內(nèi)動(dòng)態(tài)為：

（25）

2.2.3? 求解穩(wěn)定逆ud

根據(jù)式（13）求穩(wěn)定逆ud：

（26）

其中，。

綜上近似穩(wěn)定逆算法，對(duì)于復(fù)雜的過程動(dòng)力學(xué)模型或者輸出軌跡，內(nèi)動(dòng)態(tài)都可以通過增加積分迭代次數(shù)趨于收斂值。

3? 仿真驗(yàn)證

將前述反應(yīng)器模型（2）改寫為式（5）的形式，則有：

（27）

其中，模型輸入，狀態(tài)變量，模型輸出。穩(wěn)定逆方法就是已知為期望的輸出軌跡，找出輸入軌跡ud和狀態(tài)軌跡xd。

根據(jù)式（4）假設(shè)反應(yīng)器的生產(chǎn)目標(biāo)是要求產(chǎn)物濃度按照如下時(shí)變曲線變化：

（28）

為了求解進(jìn)料流量和夾套溫度調(diào)控的操作曲線即模型的穩(wěn)定逆，根據(jù)式（8）穩(wěn)定逆輸出為：

（29）

其中，是過程的不穩(wěn)定內(nèi)動(dòng)態(tài)，已知內(nèi)動(dòng)態(tài)就可以推導(dǎo)出穩(wěn)定逆ud，內(nèi)動(dòng)態(tài)方程和積分方程如下：

（30）

（31）

對(duì)于該非線性內(nèi)動(dòng)態(tài)方程使用經(jīng)典穩(wěn)定逆很難求出解析解，Picard迭代積分也無法應(yīng)用，通過利用筆者提出近似穩(wěn)定逆算法結(jié)合Picard迭代積分可以求解得到收斂的內(nèi)動(dòng)態(tài)向量，在選取n=900個(gè)分段區(qū)間進(jìn)行模型仿真，結(jié)果如圖1所示。

由圖1可知，隨著Picard迭代次數(shù)m的不斷增加，反應(yīng)器產(chǎn)物濃度相對(duì)于生產(chǎn)濃度指標(biāo)的誤差在不斷減少，這驗(yàn)證了筆者提出的近似穩(wěn)定逆算法可以與Picard迭代法結(jié)合求解非線性穩(wěn)定逆，通過模型的穩(wěn)定逆輸出即進(jìn)料流量和夾套溫度按照給定的操作曲線變化，可以完成對(duì)應(yīng)的生產(chǎn)目標(biāo)。進(jìn)料流量u1和夾套溫度u2的最佳操作曲線即Picard迭代次數(shù)為8次的近似穩(wěn)定逆輸出結(jié)果（圖2）。

為考察分段區(qū)間數(shù)對(duì)于近似穩(wěn)定逆算法的影響，對(duì)于內(nèi)動(dòng)態(tài)向量在選取n=900個(gè)和n=90兩種分段區(qū)間數(shù)進(jìn)行仿真對(duì)比，其中Picard迭代次數(shù)m=8，仿真結(jié)果如圖3所示。

由圖3可知，計(jì)算近似穩(wěn)定逆算法（n=900）模型輸出產(chǎn)物濃度y1和反應(yīng)器溫度y2相對(duì)于生產(chǎn)要求曲線的積分誤差分別為0.613 3和15.829 2；而近似穩(wěn)定逆算法（n=90）模型輸出產(chǎn)物濃度y1和反應(yīng)器溫度y2相對(duì)于生產(chǎn)要求曲線的積分誤差I(lǐng)AE分別為0.806 4和21.954 3。由此可知，分段數(shù)n越大，近似穩(wěn)定逆的誤差越小。綜合考量誤差和計(jì)算量，可以通過近似穩(wěn)定逆算法得到一個(gè)滿足生產(chǎn)過程要求的控制變量操作曲線。

4? 結(jié)束語

針對(duì)生化過程中操作曲線的求取問題，提出一種非線性穩(wěn)定逆的近似求解算法，對(duì)非線性穩(wěn)定逆算法中求解內(nèi)動(dòng)態(tài)方程方法進(jìn)行了改進(jìn)，通過將內(nèi)動(dòng)態(tài)離散化處理，利用遞推公式推導(dǎo)預(yù)覽穩(wěn)定逆中的積分方程的解，將原本求內(nèi)動(dòng)態(tài)解析解問題轉(zhuǎn)為求數(shù)值向量解問題，這樣不僅可以減少求解穩(wěn)定逆的計(jì)算量，還可以與Picard迭代法結(jié)合實(shí)現(xiàn)非線性生化過程相應(yīng)操作曲線的求解。最后通過仿真驗(yàn)證了所提算法的有效性。

需要說明的是，文中近似穩(wěn)定逆誤差大小與內(nèi)動(dòng)態(tài)跟蹤時(shí)域內(nèi)的分段數(shù)n有關(guān)，如何根據(jù)不同跟蹤目標(biāo)權(quán)衡跟蹤精度和計(jì)算量，從而選取合適的分段數(shù)n是之后需要解決的問題。

參? 考? 文? 獻(xiàn)

[1]仇力，欒小麗，劉飛.基于正交余弦相似度的間歇過程操作曲線優(yōu)化[J].化學(xué)工程，2019，47（1）：65-69.

[2]FONSECA B G，SOLEDAD M，INES C R，et al.Bioconversion in batch bioreactor of olive-tree pruning biomass optimizing treatments for ethanol production[J].Biochemical Engineering Journal，2020，164（1）：107793.

[3]THANAPIMMETHA A，SAISRIYOOT M，KHOMLAEM C，et al.A comparison of methods of ethanol production from sweet sorghum bagasse[J].Biochemical Engineering Journal，2019，151：107352.

[4]ANWESHA M，PURVI Z，RUSHIKA P，et al.Optimization of microbial fuel cell process using a novel consortium for aromatic hydrocarbon bioremediation and bioelectricity generation[J].Journal of Environmental Management，2021，298：113546.

[5]劉琳.精煉爐最佳吹氬曲線的制定與優(yōu)化[D].沈陽(yáng)：東北大學(xué)，2013.

[6]HUNT L R，MEYER G，SU R.Noncausal inverses for linear systems[J].IEEE Transactions on Automatic Control，1996，41（4）：608-611.

[7]ZOU Q Z，DEVASIA S.Preview-based stable-inversion for output tracking of linear systems[J].Journal of Dynamic Systems，Measurement，and Control，1999，121（4）：625-630.

[8]ZOU Q Z，DEVASIA S.Precision preview-based stable-inversion for nonlinear nonminimum-phase systems：The VTOL example[J].Automatica，2007，43（1）：117-127.

[9]ZOU Q Z.Optimal preview-based stable-inversion for output tracking of nonminimum-phase linear systems[J].Automatica，2009，45（1）：230-237.

[10]張有陵，劉山.基于最優(yōu)狀態(tài)轉(zhuǎn)移的非因果穩(wěn)定逆[J].控制理論與應(yīng)用，2016，33（4）：428-436.

[11]WANG H M，KIM K-S，ZOU Q Z.B-spline-decomposition-based output tracking with preview for nonminimum-phase linear systems[J].Automatica，2013，49（5）：1295-1303.

[12]ROMAGNOLI R，GARONE E.A general framework for approximated model stable inversion[J].Automatica，2019，101：182-189.

[13]NIEMIEC M P，KRAVARIS C.Nonlinear model-state feedback control for nonminimum-phase processes[J].Automatica，2003，39（7）：1295-1302.

（收稿日期：2022-10-20，修回日期：2023-07-13）