孫京誥 ，占子赫

1.華東理工大學(xué) 信息科學(xué)與工程學(xué)院，上海 200237

2.華東理工大學(xué) 化工過程先進控制和優(yōu)化技術(shù)教育部重點實驗室，上海 200237

1 引言

半間歇聚合反應(yīng)的溫度控制一直是過程控制研究的熱點問題[1-2]。非線性模型預(yù)測控制（Nonlinear Model Predictive Control，NMPC）具有顯示處理系統(tǒng)約束和在不確定環(huán)境下優(yōu)化控制的共性機理，使其成為繼PID控制之后在控制工程中被廣泛應(yīng)用和研究的先進控制技術(shù)之一[3]。NMPC運用更加合理的目標函數(shù)，比如最小化批次時間，使得反應(yīng)過程更加注重經(jīng)濟效益而不是過度地追求與目標溫度的溫差最小化[4]。然而在實際的化工過程中，往往存在實際參數(shù)與模型不匹配的情況，這使得NMPC的性能極大地下降。一類基于不確定預(yù)測模型魯棒模型預(yù)測控制旨在彌補傳統(tǒng)NMPC的不足，近年被眾多學(xué)者廣泛研究。最早有學(xué)者采用開環(huán)的min-max策略，即在通過最小化最差情況下的目標函數(shù)得到控制輸入[5]。這種方法被證實不僅具有一定程度的保守性，而且由于基于開環(huán)控制的假設(shè)與真實的閉環(huán)系統(tǒng)不匹配，導(dǎo)致其并不具備魯棒性。為了攻克這個問題，有學(xué)者提出了閉環(huán)min-max策略[6]，但是由于內(nèi)在的計算復(fù)雜度，該算法在實際應(yīng)用中的成功率有限。Tube-based mpc運用副控制器保證在不確定因素存在的情況下滿足限制條件并保證控制的穩(wěn)定性，但是該方法不能保證是最優(yōu)控制[7]。近年來，多階非線性模型預(yù)測控制（Multi-stage Nonlinear Model Predict Control，Multistage NMPC）在基于閉環(huán)魯棒MPC和運用場景樹對不確定量的建模的思想上產(chǎn)生[8-9]。在這個策略中，不確定量表示成場景樹的形式，未來的控制變量作為輔助變量來優(yōu)化，以此來降低該策略的保守性，Multi-stage NMPC已經(jīng)被證實具備很好的解決傳統(tǒng)NMPC面臨的模型參數(shù)常量與實際不匹配問題的能力。

NMPC是一類基于在線優(yōu)化的控制算法，其在線優(yōu)化計算速度直接關(guān)系到能否成功應(yīng)用。相比于線性模型預(yù)測控制，NMPC計算的最大困難在于要在有限的采樣時間內(nèi)計算一個非凸的非線性規(guī)劃問題（Nonlinear Programming，NLP）[3]。求解該問題的計算量隨決策變量的維數(shù)呈指數(shù)增長。故對于規(guī)模更大的Multi-stage NMPC，快速高效的求解方法是其能否成功應(yīng)用的關(guān)鍵。現(xiàn)有的NMPC數(shù)值算法大致分為序貫法和聯(lián)立法兩種[10]，聯(lián)立法也稱為直接轉(zhuǎn)化法，它通過同時離散狀態(tài)變量和控制變量將問題轉(zhuǎn)化為NLP。這樣將DAE或ODE系統(tǒng)的求解和優(yōu)化問題結(jié)合在一起，避免了中間解的尋找，并且具有更靈活的處理路徑不等式的能力和處理不穩(wěn)定和病態(tài)條件問題的能力[11]。有限元正交配置（Orthogonal Collocation on Finite Element，OCFE）是聯(lián)立法常用的局部離散方法，具有精度高，計算量小的特點，在動態(tài)優(yōu)化問題中有廣泛的應(yīng)用。此外，配置點的選擇合理與否，不僅關(guān)系著解的精確性，也關(guān)系著運算的復(fù)雜度，正交配置中常用的配置點有radau點、legendre點和lobatto點[12]。

本文提出了多階非線性模型預(yù)測半間歇聚合反應(yīng)控制器設(shè)計方法，將NMPC滾動優(yōu)化問題描述成最優(yōu)控制問題，其目標函數(shù)綜合考慮了溫度偏移和批次反應(yīng)時間。為了提高動態(tài)優(yōu)化問題求解效率，采用基于有限元正交配置的聯(lián)立法求解，通過實驗選取合適的radau配置點，將狀態(tài)變量和控制變量同時離散化，推導(dǎo)出基于multi-stage NMPC的可求解 NLP問題，然后借助IPOPT工具箱調(diào)用內(nèi)點法求解[13]。

2 Multi-stage NMPC

Multi-stage NMPC最重要的理論基礎(chǔ)是不確定量可以使用離散場景樹的形式表示出來，如圖1所示，在每一個節(jié)點的每一個分支都代表了一種未知的不確定因素和選擇的控制變量對系統(tǒng)的影響，這些不確定因素可能是參數(shù)誤差或者外部擾動。通過這個方程，在未來采樣時刻將得到新的信息，因此最優(yōu)化的控制變量可以當作抵消各種不確定參數(shù)的輔助變量，這種方法降低該策略的保守性。盡管在每個階段場景樹都有不同分支，但是這不表示在每個采樣點不確定量是變化的。這只是表示如果不確定量如果在第K階段是未知的，那么在第K+1階段也是未知的。因此在下一個采樣點，新的場景樹將被生成以求出下一個時刻的控制量。

圖1 Multi-stage NMPC的場景樹示意圖

由于場景樹的離散設(shè)定，在場景樹的設(shè)定下一個帶有不確定參數(shù)的不確定非線性系統(tǒng)表述如式（1）：

上式表明k時刻的狀態(tài)取決于前一個狀態(tài)，相應(yīng)的控制輸入以及不確定參數(shù)簡潔起見，所有的不確定參數(shù)和場景樹都設(shè)置為統(tǒng)一的形式，在所有的節(jié)點上都有相同數(shù)目的分支。本文約定在一個節(jié)點處，一個不確定量用三個分支來表示，假設(shè)有m個不確定量，易知在每個節(jié)點需要有3m個分支，這樣隨著預(yù)測時域的增大，優(yōu)化問題的規(guī)模也級數(shù)倍增長。為了避免場景樹的規(guī)模過大而帶來不必要的運算，設(shè)置魯棒時域NR，在魯棒時域后場景樹將不再分支，如圖1中，魯棒時域設(shè)置為2。場景樹中出現(xiàn)的(j,k)統(tǒng)一用符號S表示。每條從節(jié)點xk到葉子節(jié)點的路徑稱為一個場景。

在上述符號表述的基礎(chǔ)上，基于場景樹的最優(yōu)化問題可如式（2）描述：

公式中，N表示場景的數(shù)目，Ji和wi分別表示每個場景Si的目標方程和其相應(yīng)權(quán)重。X、U分別表示狀態(tài)變量和控制變量的取值范圍，xsoft和usoft分別表示對狀態(tài)變量和控制變量的軟化約束（Soft constraints）。式（6）稱為非預(yù)期約束條件，表示在相同的分支處，設(shè)定控制量相同。對于一般的間歇反應(yīng)模型，單個場景的目標方程有如下形式：

其中NP表示預(yù)測時域，Q和R為對角加權(quán)矩陣，目標方程的第一項表示狀態(tài)變量跟蹤設(shè)定值的懲罰項。第二項表示對控制變量變化量懲罰項，旨在防止控制變量劇烈變化。此外，限制條件包括對狀態(tài)變量和控制變量取值范圍的限制，同時，動態(tài)模型也作為最優(yōu)化問題的約束條件。

3 將動態(tài)優(yōu)化問題轉(zhuǎn)化成NLP問題

使用聯(lián)立法求解動態(tài)優(yōu)化問題需要精確的離散方法以保證對連續(xù)的軌跡的逼近精度。為了使得到的NLP問題規(guī)模不至于過大而增加求解難度，一般情況下不采用全局的離散化方法。OCFE是一種局部的離散化方法，屬于隱式Runge-Kutta法，可以提供高階的精度和良好的穩(wěn)定性。在一段采樣區(qū)間內(nèi)，定義有限元個數(shù)為nκ，有限元長度hκ。在每個有限元上通過Lagrange插值對狀態(tài)變量，控制變量進行逼近：

其中(t)表示在場景樹的(j,k)位置上，處于第κ個有限元區(qū)間內(nèi)的狀態(tài)變量值，(t)亦然。 τ∈[0,1]表示配置點。設(shè)nθ為插值的階次，則一個有限元內(nèi)的配置點 用 τ0,…,τnθ表示，本文選用Radau配置點，具體配置點信息將在第 5章給出。 lθ(τ)和 lˉθ(τ)分別表示狀態(tài)變量和控制變量的Lagrange插值多項式，可用如下形式表示：

Lagrange多項式插值具有如下性質(zhì)：

定義tk,κ,θ為第k階段第κ個有限元第θ個插值點的時間，為第k階段第κ個有限元的終端時間，該時間僅僅和有限元的個數(shù)nκ和插值階次nθ相關(guān)。于是，在場景樹中每個正交配置點上，滿足

于是可以得到離散化后的NLP問題：

其中，式（12）通過將每個有限元微分得到，其中微分變量需要保持狀態(tài)的連續(xù)性，式（13）保證了每個相鄰有限元之間的連續(xù)性，式（14）保證了相鄰控制時域內(nèi)的連續(xù)性。狀態(tài)變量約束和控制變量約束條件由式（15）給出。

4 半間歇聚合反應(yīng)模型

為了驗證提出的算法的可行性，本文對一種間歇性化工聚合反應(yīng)進行研究，該模型來自于一個實際的化工間歇聚合反應(yīng)模型，它的機理模型由BASF SE公司提供。

4.1 機理模型

該系統(tǒng)主要由一個反應(yīng)釜，夾套以及一個外部熱轉(zhuǎn)換器（External Heat Exchanger，EHE）組成，夾套和EHE可以用來控制反應(yīng)釜內(nèi)的溫度。其工藝圖如圖2。在一個批次的反應(yīng)中，開始階段將反應(yīng)所需的引發(fā)劑，水等反應(yīng)所需物質(zhì)加入反應(yīng)釜內(nèi)，將釜溫加熱到設(shè)定溫度后進入進料階段，即開始持續(xù)向釜內(nèi)加入單體，當加入的單體達到一定質(zhì)量時停止進料，進入保持階段，直到產(chǎn)生到規(guī)定質(zhì)量的聚合物后反應(yīng)停止。

圖2 半間歇式聚合反應(yīng)工藝圖

根據(jù)物料守恒和能量守恒，其機理模型由主要由8個常微分方程構(gòu)成，如式（17）～（24）所示：

該機理模型包括水mW，單體mA和聚合產(chǎn)物mP的物料守恒方程以及反應(yīng)釜內(nèi)溫TR，管道溫度TS，夾套溫度TM，外部熱轉(zhuǎn)換器溫度TEK和冷卻劑在熱轉(zhuǎn)換器的出口溫度TAWT的能量守恒方程。反應(yīng)過程中有三個控制變量，分別是進料速率mF，外部熱交換器入口溫度夾套入口冷卻劑溫度TINM。

U表示反應(yīng)釜內(nèi)聚合物與單體的比例：

mges表示釜內(nèi)物質(zhì)的總質(zhì)量：

kR1、kR2分別表示表示釜內(nèi)的反應(yīng)速率和外部熱轉(zhuǎn)換器的反應(yīng)速率，如式（27）、（28）所示：

kK表示反應(yīng)釜內(nèi)的混合物總傳熱系數(shù)，由式（29）給出：

mA,R表示釜內(nèi)的單體質(zhì)量，由式（30）得出：

4.2 控制目標

在一個批次反應(yīng)過程中，在滿足安全條件限制和產(chǎn)物質(zhì)量的前提下，要在盡可能短的時間內(nèi)通過調(diào)整mF、TINAWT、TINM使反應(yīng)安全穩(wěn)定迅速地進行。溫度的控制在聚合反應(yīng)中占據(jù)了很重要的位置，它和聚合產(chǎn)物的質(zhì)量息息相關(guān)。在本問題中，反應(yīng)釜的溫度應(yīng)該被控制在設(shè)定溫度Tset=(363.15±2.0)K范圍內(nèi)。在真實的反應(yīng)過程中，出于對安全性和反應(yīng)器材耐受能力的考慮，反應(yīng)過程溫度還需要滿足安全限制。在本模型中，為了避免冷卻失敗的情況，反應(yīng)釜最高溫度被限制在382.15 K。在任意時刻反應(yīng)釜的安全溫度由式（31）給出：

此外，加入反應(yīng)釜的單體的最大質(zhì)量為30 000 kg。進料階段投入單體，當添加完所需的所有的單體，反應(yīng)進入保持階段，如果聚合物的質(zhì)量達到20 680 kg，則反應(yīng)停止。若在給定的單體范圍內(nèi)沒有生產(chǎn)出規(guī)定質(zhì)量的聚合物則反應(yīng)失敗。在真實的系統(tǒng)中，通常模型的參數(shù)不能被準確地測量。在本模型中，有兩個關(guān)鍵的參數(shù)：反應(yīng)熱焓ΔHR和反應(yīng)速率k0，這兩個參數(shù)通常會在其給定值的±30%內(nèi)變化，且在一個批次的反應(yīng)過程中值不變。模型的狀態(tài)變量的初始值與取值范圍和控制變量的取值范圍分別如表1，表2所示。

表1 狀態(tài)變量初始值和取值范圍

表2 控制變量取值范圍

5 仿真分析

所有的計算都運行在Windows環(huán)境臺式電腦（Inter Core i7-6700K@4.00 GHz 16 GB DDR4 2 801 MHz）。Multi-stage NMPC所產(chǎn)生的NLP問題都采用最優(yōu)化工具包IPOPT來求解，IPOPT工具包內(nèi)部使用內(nèi)點法對最優(yōu)化問題進行求解。設(shè)定Multi-stage NMPC和NMPC的采樣時間tstep=50 s，預(yù)測時域Np=10，設(shè)定Multi-stage NMPC的魯棒時域NR=1。反應(yīng)模型中兩個不確定反應(yīng)熱焓ΔHR和反應(yīng)速率k0在其給定值的±30%內(nèi)變化，其真實值分表由，表示。由于有兩個不確定量，在魯棒時域內(nèi)，場景樹的每個節(jié)點都有9個分支，分支的不確定參數(shù)以最壞情況和標準情況的組合，即=(0.7,1.0,1.3)×ΔHR分別與=(0.7,1.0,1.3)×k0組合得到的9種可能。優(yōu)化問題目標方程綜合考慮了溫度跟蹤和反應(yīng)的速率，設(shè)定式（6）中-1,100,0,0,0,0,0,0],R=diag[0.2,0.03,0.02]。此外，仿真過程中假定所需要的變量都可以被準確地測量。

5.1 正交配置點的選取

配置點τ的選擇對求解的精度，速度以及解的存在與否都有很大的影響[15]。若采用等間隔方式選取配置點，當點過于密集時可能產(chǎn)生龍格（Runge）現(xiàn)象[14]。本文選擇radau點作為配置點，使用radau配點的OCFE可以有效避免龍格現(xiàn)象，并獲得平滑的逼近軌跡，且相比于legendre點，radau點可以提供更佳的穩(wěn)定性。radau配置點取值如表3。

表3 radau配置點取值

設(shè)定有限元個數(shù)nκ=2，配置點階次nθ分別取2、3、4進行實驗仿真，由于模型參數(shù)的不確定性，在ΔHrealR=(0.7,1.0,1.3)×ΔHR分別與 kreal0=(0.7,1.0,1.3)×k0組合得到的9種可能情況下進行實驗，分別得到批次反應(yīng)時間和求解一次優(yōu)化問題的平均時間，并以箱線圖的形式畫出，如圖3、4。從圖3中可以看出，隨著nθ的增加，批次反應(yīng)時間有減少的趨勢，這符合目標方程的設(shè)定?？梢妌θ在一定范圍內(nèi)增加可以通過更好地逼近真實曲線而帶來更好的控制效果。而從圖4反應(yīng)出，隨著nθ的增加，在反應(yīng)過程中，在每個采樣時間求解最優(yōu)化問題所需的平均運算時間也大幅度地增加，以中位數(shù)為例，nθ=2時，中位數(shù)為2.1 s，而在 nθ=3,4 時，這一數(shù)值分別達到了6.5 s和26.2 s。故雖然配點階次增加可以帶來控制效果的略微提升，但由于過度地增加了優(yōu)化問題的計算時間使得不符合實際過程控制的實時性要求。

圖3 radau配置點階次為2、3、4時，一個批次反應(yīng)所需時間

圖4 radau配置點階次為2、3、4時，求解一次優(yōu)化問題平均時間

5.2 控制效果分析

根據(jù)5.1節(jié)實驗結(jié)果，設(shè)定nκ=2，nθ=2，分別采用直接radau配置的Multi-stage NMPC（以下簡稱為RMSNMPC）和傳統(tǒng)NMPC對模型進行仿真，模型參數(shù)設(shè)置如5.1節(jié)。在9種參數(shù)組合中有3種情況NMPC無法完成一個批次的反應(yīng)，故不作圖分析。其他仿真結(jié)果如圖5～9，圖中的紅線表示溫限。從圖中可以看出，在任意一個情形下，RMSNMPC均可以較好地完成一個批次的反應(yīng)，釜溫不僅被控制在安全溫度以內(nèi)，也被控制在(363.15±2.0)K精度范圍內(nèi)。而NMPC在多數(shù)情形下，如圖6、7、9，反應(yīng)釜釜溫都沒有達到要求的控制精度要求。更為嚴重的是，在 ΔHrealR=1.0×ΔHR，kreal0=1.3×k0的情況下，如圖7，反應(yīng)進行到60 min時，處于進料階段，釜溫一度超過了安全溫度，這會給實際生產(chǎn)過程帶來嚴重的安全隱患。此外，雖然NMPC在標稱情況和 ΔHrealR=0.7×ΔHR，kreal0=1.3×k0的情況下可以達到溫度的精度和安全要求，但是進料量都有劇烈震蕩的現(xiàn)象發(fā)生，這會給實際控制帶來操作難度。相比之下，RMSNMPC在各個情形下都可以帶來更加平滑的控制量?？梢娤啾扔趥鹘y(tǒng)的NMPC，RMSNMPC利用場景樹模擬不確定量的優(yōu)勢，可以很好地處理具有不確定參數(shù)的模型。

圖5 =1.0×ΔHR，=1.0×k0時，釜溫、安全溫度、進料量、夾套入口冷卻劑溫度曲線圖

圖6 =1.0×ΔHR，=1.3×k0時，釜溫、安全溫度、進料量、夾套入口冷卻劑溫度曲線圖

圖7 =1.3×ΔHR，=1.0×k0時，釜溫、安全溫度、進料量、夾套入口冷卻劑溫度曲線圖

圖8 =0.7×ΔHR，=1.3×k0時，釜溫、安全溫度、進料量、夾套入口冷卻劑溫度曲線圖

圖9 =1.3×ΔHR=0.7×k0時，釜溫、安全溫度、進料量、夾套入口冷卻劑溫度曲線圖

6 結(jié)束語

針對半間歇聚合反應(yīng)模型中存在不確定參數(shù)的情況，基于場景樹模擬不確定量的Multi-stage NMPC方法，提出使用基于有限元正交配置的聯(lián)立法對優(yōu)化問題進行求解。通過大量實驗選取最佳的radau配置點，形成有效的RMSNMPC方法。仿真結(jié)果表明在標稱模型和具有不確定參數(shù)的模型下，所采用的算法不僅能使反應(yīng)滿足被控的溫度精度要求和安全要求，還具備良好實時求解能力和魯棒性。相比于傳統(tǒng)的NMPC算法，更適用于復(fù)雜的半間歇聚合反應(yīng)的控制。

[1]Wang J，Wang Y，Cao L，et al.Adaptive interative learning control based on unfalsified strategy for Chylla-Haase reactor[J].IEEE/CAA Journal of Automatica Sinica，2014，1（4）：347-360.

[2]Jang H，Lee J H，Biegler L T.A robust NMPC scheme for semi-batch polymerization reactors[J].IFAC-Papers on Line，2016，49（7）：37-42.

[3]何德峰，丁寶蒼，于樹友.非線性系統(tǒng)模型預(yù)測控制若干基本特點與主題回顧[J].控制理論與應(yīng)用，2013，30（3）：273-287.

[4]Idris E A N，Engell S.Economics-based NMPC strategies for the operation and control of a continuous catalytic distillation process[J].Journal of Process Control，2012，22（10）：1832-1843.

[5]Campo P J，Morari M.Robust model predictive control[C]//American Control Conference，1987：1021-1026.

[6]Lee J H，Yu Z.Worst-case formulations of model predictive control for systems with bounded parameters[J].Automatica，1997，33（5）：763-781.

[7]Mayne D Q，Kerrigan E C，Van Wyk E J，et al.Tube based robust nonlinear model predictive control[J].International Journal of Robust and Nonlinear Control，2011，21（11）：1341-1353.

[8]Lucia S，F(xiàn)inkler T，Engell S.Multi-stage nonlinear model predictive control applied to a semi-batch polymerization reactor under uncertainty[J].Journal of Process Control，2013，23（9）：1306-1319.

[9]Lucia S，Engell S.Multi-stage and two-stage robust nonlinear model predictive control[C]//Nonlinear Model Predictive Control，2012，4（1）：181-186.

[10]陳楊，邵之江，錢積新，等.聯(lián)立法中全局和局部正交配置算法[J].化工學(xué)報，2010（2）：384-391.

[11]黃森.動態(tài)優(yōu)化問題的移動有限元算法研究[D].杭州：浙江大學(xué)，2012.

[12]Canuto C，Hussaini M Y，Quarteroni A M，et al.Spectral methods in fluid dynamics[M].New York：Springerverlay，2012.

[13]Biegler L T.An overview of simultaneous strategies for dynamic optimization[J].Chemical Engineering and Processing：Process Intensification，2007，46（11）：1043-1053.

[14]Ascher U M，Petzold L R.Computer methods for ordinary differential equations and differential-algebraic equations[M].Philadephia，PA，USA：SIAM，1998.

[15]陳楊.基于微分—代數(shù)混合方程機理模型的非線性預(yù)測控制[D].杭州：浙江大學(xué)，2011.