劉星偉，賈勝坤，羅祎青，袁希鋼，2

（1 天津大學(xué)化工學(xué)院，天津 300354； 2 化學(xué)工程聯(lián)合國(guó)家重點(diǎn)實(shí)驗(yàn)室（天津大學(xué)），天津 300354）

引 言

分離過(guò)程是化學(xué)工業(yè)生產(chǎn)過(guò)程中的重要過(guò)程，而精餾分離是化工生產(chǎn)過(guò)程中最主要的分離方法[1]。精餾塔的嚴(yán)格優(yōu)化與設(shè)計(jì)對(duì)于降低精餾過(guò)程能耗，節(jié)約生產(chǎn)成本，減少生產(chǎn)過(guò)程的碳排放和環(huán)境污染具有重要意義[2-4]。

精餾塔優(yōu)化設(shè)計(jì)的決策變量既包含回流比、壓力、再沸比等連續(xù)變量，又包含總塔板數(shù)、進(jìn)料板位置等離散變量，因此對(duì)精餾塔決策變量進(jìn)行同時(shí)優(yōu)化屬于混合整數(shù)非線性規(guī)劃(mixed-integer nonlinear programming, MINLP)問(wèn)題，其優(yōu)化求解面臨很多困難。在傳統(tǒng)的精餾塔優(yōu)化工作中，往往直接將總塔板數(shù)和進(jìn)料位置作為整數(shù)變量進(jìn)行優(yōu)化[5-7]。Viswanathan 等[8]假設(shè)回流流股和再沸流股同時(shí)進(jìn)入多個(gè)塔板，并且采用0-1二元變量來(lái)表示流股具體進(jìn)入的位置，構(gòu)造了一種精餾塔級(jí)數(shù)優(yōu)化的超結(jié)構(gòu)模型，并在許多工作中都得到了廣泛應(yīng)用[9-11]。但是，對(duì)于此模型的優(yōu)化仍然是MINLP 問(wèn)題，后續(xù)研究者提出了多種方法，將0-1二元變量進(jìn)行連續(xù)化松弛，以達(dá)到降低優(yōu)化問(wèn)題求解難度的目的。Lang等[12]假設(shè)流股可進(jìn)入多個(gè)平衡級(jí)，采用可微分的微分分布函數(shù)來(lái)表示進(jìn)入一個(gè)平衡級(jí)的流量占總流量的比例，從而使進(jìn)料位置和總塔板數(shù)等整數(shù)變量轉(zhuǎn)化為連續(xù)變量，實(shí)現(xiàn)了塔板數(shù)的連續(xù)化。Neves 等[13]受到Lang 等工作[12]的啟發(fā)，不再采用微分分布函數(shù)，而是直接采用0-1 之間的分?jǐn)?shù)代表進(jìn)入一個(gè)平衡的流量占總流股流量的比例，再通過(guò)不等式約束使分?jǐn)?shù)趨近整數(shù)，最終得到整數(shù)。Kraemer等[14]同樣首先假設(shè)回流流股和再沸流股同時(shí)進(jìn)入多個(gè)塔板，并且采用0-1二元變量來(lái)表示流股具體進(jìn)入的位置，然后再松弛二元變量為連續(xù)變量，最后引入非線性Fischer-Burmeister 函數(shù)來(lái)驅(qū)使這些連續(xù)變量收斂為0-1二元變量值。

在上述精餾塔的優(yōu)化工作中，研究者通過(guò)使進(jìn)料流股和回流流股同時(shí)進(jìn)入多個(gè)平衡級(jí)，建立進(jìn)料位置和總塔板數(shù)優(yōu)化的超結(jié)構(gòu)。這種精餾塔模型的超結(jié)構(gòu)較為復(fù)雜，因此很難應(yīng)用到復(fù)雜精餾流程中。本文從結(jié)構(gòu)優(yōu)化的角度，將精餾塔優(yōu)化中較為困難的進(jìn)料位置和總塔板數(shù)優(yōu)化轉(zhuǎn)化為空間結(jié)構(gòu)的分布優(yōu)化，即在精餾段和提餾段預(yù)設(shè)多于需要的平衡級(jí)數(shù)，采用結(jié)構(gòu)參數(shù)代表平衡級(jí)是否存在，再通過(guò)精餾段和提餾段的實(shí)際平衡級(jí)數(shù)來(lái)表示進(jìn)料位置和總塔板數(shù)。

在結(jié)構(gòu)優(yōu)化的研究中，信賴域優(yōu)化算法更常被用來(lái)求解優(yōu)化問(wèn)題。不同于線搜索優(yōu)化算法，信賴域優(yōu)化算法采用先步長(zhǎng)、后方向的策略，其在當(dāng)前迭代點(diǎn)的鄰域內(nèi)構(gòu)造優(yōu)化子問(wèn)題，并逐步逼近原問(wèn)題的最優(yōu)解。信賴域算法具有很好的適用性和穩(wěn)健性，并且可以保證優(yōu)化的全局收斂性[15]，另外，信賴域算法可以方便地與流程模擬軟件進(jìn)行聯(lián)用，因此，其在精餾塔優(yōu)化中有許多優(yōu)勢(shì)。

本文從結(jié)構(gòu)優(yōu)化的角度對(duì)精餾塔進(jìn)行建模，采用信賴域優(yōu)化算法對(duì)模型進(jìn)行優(yōu)化，并搭建了流程模擬軟件Aspen Custom Modeler 與Matlab 結(jié)合的優(yōu)化方法。為了證明提出的優(yōu)化方法的可行性和收斂性，本文采用該優(yōu)化方法分別對(duì)單一精餾塔、直接精餾序列和部分熱耦合精餾塔進(jìn)行優(yōu)化設(shè)計(jì)。

1 模型與優(yōu)化算法

1.1 精餾模型

本文從結(jié)構(gòu)優(yōu)化的角度建立精餾模型，將總塔板數(shù)和進(jìn)料位置的優(yōu)化轉(zhuǎn)化為塔板空間分布的結(jié)構(gòu)優(yōu)化，并用0-1二元結(jié)構(gòu)參數(shù)表示塔板的存在性。如圖1 所示，在精餾段和提餾段預(yù)先設(shè)置多于需要的平衡級(jí)數(shù)，每一個(gè)方框代表一個(gè)平衡級(jí)，當(dāng)二元結(jié)構(gòu)參數(shù)為1時(shí)，代表平衡級(jí)起到分離作用（灰色），取0 則表示平衡級(jí)不存在。這樣，通過(guò)優(yōu)化0-1 二元結(jié)構(gòu)參數(shù)的分布，就可以實(shí)現(xiàn)對(duì)總塔板數(shù)和進(jìn)料位置的同時(shí)優(yōu)化，將精餾過(guò)程的優(yōu)化轉(zhuǎn)化為二元結(jié)構(gòu)參數(shù)分布結(jié)構(gòu)的優(yōu)化。

圖1 精餾塔結(jié)構(gòu)參數(shù)模型示意圖Fig.1 Schematic diagram of distillation column structure parameter model

為了避免優(yōu)化求解0-1 二元變量的MINLP 問(wèn)題，本文采用Dowling 等[16]提出的繞流效率模型將二元變量連續(xù)化。如圖2 所示，圖中L為液相流股，V為氣相流股，其基本思想是為每一平衡級(jí)增加繞流流股，流入第j個(gè)平衡級(jí)的氣液相流股，被分為兩部分：一部分進(jìn)入該平衡級(jí)進(jìn)行分離；另一部分繞過(guò)該平衡級(jí)。進(jìn)入平衡級(jí)的流量占總流量的比例則用繞流效率εj表示。再將該平衡級(jí)上的汽液平衡流股與繞流流股混合，就得到離開該平衡級(jí)的氣液相流股。

圖2 繞流效率模型Fig.2 Illustration of bypass efficiency model

1.2 移動(dòng)漸近線算法

在優(yōu)化算法方面，本文采用信賴域優(yōu)化算法中的移動(dòng)漸近線法（MMA），其最早由Svanberg[17-18]提出，并被應(yīng)用于機(jī)械工程的結(jié)構(gòu)優(yōu)化問(wèn)題中[19-20]。該算法對(duì)于一個(gè)一般的非線性約束優(yōu)化問(wèn)題[式(1)]，在當(dāng)前迭代點(diǎn)xk處采用子問(wèn)題[式(2)]進(jìn)行近似。

圖3 MMA算法Fig.3 MMA algorithm

1.3 優(yōu)化方法

MMA 算法僅需要目標(biāo)函數(shù)值、約束函數(shù)值、目標(biāo)函數(shù)和約束函數(shù)關(guān)于決策變量的梯度信息即可構(gòu)造出子問(wèn)題，求解子問(wèn)題即可得到下一個(gè)迭代點(diǎn)，因此很容易和流程模擬軟件結(jié)合使用。本文在Aspen Custom Modeler（ACM）中建立精餾模型方程，在Matlab 中執(zhí)行MMA 優(yōu)化算法，通過(guò)ACM 和Matlab 聯(lián)用傳輸模型和優(yōu)化算法之間的信息，具體的優(yōu)化方法如圖4所示。

圖4 精餾塔優(yōu)化方法Fig.4 Optimization method of distillation column

本文采用聯(lián)立方程法進(jìn)行精餾模型的求解，該方法可以較容易地獲取優(yōu)化所需要的梯度信息。但是，由于精餾塔的MESH 方程是一個(gè)強(qiáng)耦合的非線性方程組，為了解決其穩(wěn)態(tài)解的求解困難問(wèn)題，本文采用Ma 等[21]提出的虛擬瞬態(tài)連續(xù)性模型(PTC)來(lái)輔助精餾模型方程的求解。該模型引入虛擬持液量的概念，采用動(dòng)態(tài)模擬方程來(lái)實(shí)現(xiàn)MESH 方程組的求解。當(dāng)該動(dòng)態(tài)模型方程積分至穩(wěn)態(tài)時(shí)，即可以獲得MESH 方程的穩(wěn)態(tài)解。ACM 可以方便地輸出優(yōu)化所需要的信息，并將優(yōu)化信息傳入Matlab 執(zhí)行MMA 的一步迭代過(guò)程。之后，優(yōu)化算法將會(huì)給出下一個(gè)迭代點(diǎn)，將此迭代點(diǎn)傳入ACM 中進(jìn)行穩(wěn)態(tài)模擬計(jì)算。若穩(wěn)態(tài)模擬計(jì)算失敗，則采用虛擬瞬態(tài)連續(xù)性方程進(jìn)行動(dòng)態(tài)積分，進(jìn)而獲得當(dāng)前迭代點(diǎn)的穩(wěn)態(tài)模擬結(jié)果。因?yàn)閯?dòng)態(tài)積分過(guò)程所需時(shí)間較長(zhǎng)，因此優(yōu)先執(zhí)行穩(wěn)態(tài)模擬過(guò)程，穩(wěn)態(tài)模擬失敗再執(zhí)行動(dòng)態(tài)積分過(guò)程，以減少模擬所需要的時(shí)間。重復(fù)以上迭代過(guò)程，直到滿足收斂條件便可獲得最終的優(yōu)化結(jié)果。

2 算例分析

本文采用1.3 節(jié)中所提出的優(yōu)化方法對(duì)3 個(gè)算例進(jìn)行優(yōu)化，優(yōu)化的目標(biāo)函數(shù)均為年度總費(fèi)用（TAC）。TAC 由設(shè)備費(fèi)和操作費(fèi)構(gòu)成，設(shè)備費(fèi)根據(jù)精餾塔塔板數(shù)、塔徑和換熱器的換熱面積計(jì)算，操作費(fèi)由再沸器和冷凝器的換熱量算出，具體計(jì)算公式可參考Jiang 等[22]的工作。每一個(gè)優(yōu)化算例分別從5 個(gè)不同的初始點(diǎn)進(jìn)行優(yōu)化，以此說(shuō)明所提優(yōu)化方法的可行性及收斂性。

2.1 乙苯-苯乙烯二元混合分離

在本算例中，一個(gè)單精餾塔被用于分離乙苯和苯乙烯的混合物，混合物的進(jìn)料情況如圖5 所示，規(guī)定塔頂和塔釜產(chǎn)品的摩爾分?jǐn)?shù)高于95%。因?yàn)楸揭蚁┰趬毫^高時(shí)易于發(fā)生聚合，因此塔頂?shù)膲毫x為6 kPa。優(yōu)化變量包括每一平衡級(jí)的繞流效率、回流比（RR）和再沸器的汽化分率（Vf）。本算例為每個(gè)塔節(jié)設(shè)置25 個(gè)平衡級(jí)，每個(gè)塔節(jié)有25 個(gè)繞流效率作為決策變量，因此該優(yōu)化問(wèn)題的決策變量數(shù)為52，約束個(gè)數(shù)為2。結(jié)合目標(biāo)函數(shù)和精餾過(guò)程應(yīng)滿足的約束條件，可以得到式（3）所示的約束問(wèn)題，該問(wèn)題的結(jié)構(gòu)同樣適用于后面兩個(gè)優(yōu)化算例。

圖5 乙苯-苯乙烯分離塔的結(jié)構(gòu)Fig.5 Structure of ethylbenzene-styrene column

本文的優(yōu)化初始值選擇回流比為5，再沸器的汽化分率為0.5，繞流效率分別從5 個(gè)不同的初始點(diǎn)0.1、0.3、0.5、0.7 和0.9 進(jìn)行優(yōu)化，不同初始點(diǎn)的優(yōu)化結(jié)果如表1 所示。表中的N1和N2分別表示精餾段和提餾段繞流效率之和，即精餾段和提餾段的塔板數(shù)。

表1 乙苯-苯乙烯分離塔的優(yōu)化結(jié)果Table 1 Optimization results of ethylbenzene-styrene column

如表1所示，在5個(gè)不同繞流效率初始點(diǎn)的優(yōu)化計(jì)算中，繞流效率ε都收斂到0 或1，因此，最終得到的塔板數(shù)為整數(shù)。繞流效率模型在不對(duì)分?jǐn)?shù)繞流效率進(jìn)行懲罰時(shí)，繞流效率仍趨于0 或1，Dowling 等[16]解釋這一現(xiàn)象，在繞流效率為分?jǐn)?shù)時(shí)，意味著存在一部分流股繞過(guò)塔板，此時(shí)的熱力學(xué)效率較低，因此繞流效率趨于0 或1。許多應(yīng)用繞流效率模型對(duì)精餾塔進(jìn)行優(yōu)化的工作中，最終都得到了整數(shù)繞流效率的優(yōu)化結(jié)果[21,23-24]。該結(jié)果一方面說(shuō)明了繞流效率模型的有效性；另一方面也說(shuō)明移動(dòng)漸近線法可有效找到該問(wèn)題的最優(yōu)值。另外，從不同的初始值進(jìn)行優(yōu)化，MMA算法均可找到該問(wèn)題的一個(gè)局部最優(yōu)值，說(shuō)明MMA 優(yōu)化算法的收斂性較好。圖6給出了以繞流效率0.1 為初值時(shí)，N1（精餾段繞流效率之和）、N2（提餾段繞流效率之和）和目標(biāo)函數(shù)TAC隨迭代次數(shù)的變化。從圖中可以看出，僅需要36次迭代即可得到最終的優(yōu)化結(jié)果。

圖6 N1、N2和TAC隨迭代次數(shù)的變化Fig.6 Variation of N1,N2 and TAC with iteration

從5 個(gè)不同的初始點(diǎn)進(jìn)行優(yōu)化，最終優(yōu)化結(jié)果并不完全相同，說(shuō)明繞流效率模型有較多的局部最優(yōu)解。通過(guò)枚舉可以知道，該精餾結(jié)構(gòu)的最優(yōu)精餾段平衡級(jí)數(shù)為19，提餾段塔板數(shù)為18，回流比為6.48，塔釜汽化分率為0.8380，對(duì)應(yīng)的最優(yōu)TAC 為324.474×104USD/a。對(duì)比得到的5 個(gè)最優(yōu)結(jié)果，雖然均未得到全局最優(yōu)值，但與全局最優(yōu)值相差較小，這說(shuō)明MMA 優(yōu)化算法用于精餾塔優(yōu)化同樣可以獲得較好的結(jié)果。本方法無(wú)法保證獲得全局最優(yōu)解，但是具有較好的穩(wěn)健性，為獲得較好的局部最優(yōu)解甚至全局最優(yōu)解，可以從多個(gè)初始點(diǎn)出發(fā)進(jìn)行計(jì)算，以選取最好的局部最優(yōu)解。

2.2 直接精餾塔序列

直接精餾塔序列是分離三組元混合物的基本方法之一。在本算例中，采用直接序列對(duì)苯、甲苯和對(duì)二甲苯進(jìn)行分離。產(chǎn)品的進(jìn)料情況如圖7 所示。本算例要求三個(gè)產(chǎn)品的摩爾分?jǐn)?shù)不低于95%，塔頂?shù)膲毫潭? bar（1 bar=105Pa），優(yōu)化變量包括每個(gè)平衡級(jí)的繞流效率、兩個(gè)精餾塔的回流比（RR1、RR2）和兩個(gè)精餾塔再沸器汽化分率（Vf1、Vf2）。本算例仍為每個(gè)塔節(jié)設(shè)置25個(gè)平衡級(jí)，因此本優(yōu)化算例包括3 個(gè)約束條件，104 個(gè)優(yōu)化變量。表2 展示了最終的優(yōu)化結(jié)果，表中的N1、N2、N3和N4的含義如圖7所示。

圖7 直接序列精餾塔結(jié)構(gòu)Fig.7 Structure of direct sequence distillation column

如表2 所示，在5 個(gè)不同初始點(diǎn)的優(yōu)化計(jì)算中，從不同初始值開始優(yōu)化，繞流效率最終均收斂為0或1，相應(yīng)的平衡級(jí)數(shù)為整數(shù)，這說(shuō)明在本算例中MMA 優(yōu)化算法同樣可以找到局部最優(yōu)解。本算例優(yōu)化變量比2.1 節(jié)算例多一倍，在這種情況下，不同初始值均能收斂到局部最優(yōu)解，說(shuō)明MMA 優(yōu)化算法的收斂性和穩(wěn)定性較好。

表2 直接序列精餾塔的優(yōu)化結(jié)果Table 2 Optimization results of direct sequence distillation column

2.3 部分熱耦合的精餾過(guò)程

精餾操作所消耗的能量占化工生產(chǎn)總能量消耗的很大比例，部分熱耦合的精餾過(guò)程常常被作為一種精餾設(shè)計(jì)的節(jié)能措施[25-30]。但是，由于部分熱耦合精餾流程中存在質(zhì)量和能量的耦合，因此求解難度增加。本節(jié)考慮對(duì)圖8 的精餾結(jié)構(gòu)進(jìn)行優(yōu)化，進(jìn)料為正己烷、正庚烷和正辛烷。兩個(gè)精餾塔的塔頂壓力均固定為1 bar，三個(gè)產(chǎn)品純度均要求摩爾分?jǐn)?shù)不低于95%。優(yōu)化變量包括每個(gè)平衡級(jí)的繞流效率、兩個(gè)精餾塔的回流比（RR1、RR2）、第二個(gè)精餾塔的再沸器汽化分率（Vf）以及進(jìn)入第二個(gè)精餾塔精餾段的氣量占第二個(gè)精餾塔提餾段氣量的比例（Sf）。本算例仍為每個(gè)塔節(jié)設(shè)置25 個(gè)平衡級(jí)，因此本優(yōu)化算例包括3 個(gè)約束條件，104 個(gè)優(yōu)化變量。表3 展示了最終的優(yōu)化結(jié)果，表中的N1、N2、N3和N4的含義如圖8所示。

圖8 部分熱耦合精餾結(jié)構(gòu)Fig.8 Structure of partially thermally coupled distillation process

如表3 所示，在5 個(gè)不同初始點(diǎn)的優(yōu)化計(jì)算中，不同初始值的優(yōu)化結(jié)果的繞流效率均為0或1，這說(shuō)明繞流效率模型在具有熱耦合精餾結(jié)構(gòu)中依然有效。5 個(gè)不同的初始值均獲得了優(yōu)化結(jié)果，這一方面說(shuō)明所提優(yōu)化框架穩(wěn)健性較好，對(duì)于含有熱耦合結(jié)構(gòu)的優(yōu)化過(guò)程依然有效，另一方面說(shuō)明MMA 優(yōu)化算法有較好的收斂性，在含有熱耦合結(jié)構(gòu)的優(yōu)化中依然有效。

表3 部分熱耦合精餾過(guò)程優(yōu)化結(jié)果Table 3 Optimization results of partially thermally coupled distillation process

3 結(jié) 論

（1）從結(jié)構(gòu)優(yōu)化的角度，將精餾塔總塔板數(shù)和進(jìn)料位置的優(yōu)化轉(zhuǎn)化為二元結(jié)構(gòu)參數(shù)分布結(jié)構(gòu)的優(yōu)化，并結(jié)合繞流效率模型將二元離散變量轉(zhuǎn)化為連續(xù)變量。該模型可以降低精餾塔模型的復(fù)雜程度，實(shí)現(xiàn)對(duì)總塔板數(shù)和進(jìn)料位置的同時(shí)優(yōu)化。

（2）采用聯(lián)立方程法對(duì)精餾塔進(jìn)行建模，采用虛擬瞬態(tài)連續(xù)性方程來(lái)輔助精餾穩(wěn)態(tài)模擬計(jì)算，實(shí)現(xiàn)了精餾塔決策變量的同時(shí)優(yōu)化。分別對(duì)一個(gè)單一精餾塔、直接序列精餾塔和熱耦合結(jié)構(gòu)精餾塔進(jìn)行優(yōu)化設(shè)計(jì)。結(jié)果表明，虛擬瞬態(tài)連續(xù)性方程可有效輔助穩(wěn)態(tài)模擬的計(jì)算。

（3）結(jié)合MMA 算法的特點(diǎn)，將Aspen Custom Modeler 和Matlab 耦合鏈接，充分利用了Aspen Custom Modeler 的模擬能力和Matlab 的優(yōu)化計(jì)算能力。該優(yōu)化方法可保證本文3個(gè)優(yōu)化算例的正常運(yùn)行，體現(xiàn)了該方法的穩(wěn)健性。

從優(yōu)化算例來(lái)看，初始值的選取會(huì)對(duì)優(yōu)化結(jié)果產(chǎn)生影響，如何選取初始值以獲得更好的結(jié)果值得進(jìn)一步探究。本文建立的優(yōu)化方法對(duì)許多優(yōu)化算法同樣適用，因此，不同優(yōu)化算法在精餾塔優(yōu)化中的表現(xiàn)是未來(lái)一個(gè)值得研究的方向。