戴吟 劉飛

基金項目：國家自然科學基金項目（批準號：61833007）資助的課題。

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

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

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

DOI：10.20030/j.cnki.1000-3932.202306000

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

關鍵詞? 操作曲線? 動力學模型? 非線性穩(wěn)定逆? 內(nèi)動態(tài)方程? 化學反應器

中圖分類號? TP273? ? ? ? ?文獻標志碼? A? ? ? ? ?文章編號? 1000-3932（2023）06-0000-00

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

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

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

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

1? 操作曲線求取問題

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

（1）

反應器的動力學模型包括物料A、B的摩爾平衡和能量平衡：

（2）

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

（3）

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

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

（4）

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

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

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

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

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

（5）

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

設為期望的輸出軌跡，需要找出對應的輸入軌跡和狀態(tài)軌跡，即：

（6）

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

（7）

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

（8）

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

（9）

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

（10）

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

（11）

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

（12）

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

（13）

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

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

（14）

（15）

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

（16）

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

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

（17）

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

（18）

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

（19）

當?shù)螖?shù)時，。

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

（20）

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

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

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

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

（21）

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

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

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

（22）

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

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

（23）

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

（24）

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

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

（25）

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

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

（26）

其中，。

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

3? 仿真驗證

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

（27）

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

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

（28）

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

（29）

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

（30）

（31）

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

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

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

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

4? 結(jié)束語

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

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

參? 考? 文? 獻

[1]仇力，欒小麗，劉飛.基于正交余弦相似度的間歇過程操作曲線優(yōu)化[J].化學工程，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].沈陽：東北大學，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].控制理論與應用，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）