基于信賴域算法的精餾塔優化

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

（1 天津大學化工學院，天津 300354； 2 化學工程聯合國家重點實驗室（天津大學），天津 300354）

引 言

分離過程是化學工業生產過程中的重要過程，而精餾分離是化工生產過程中最主要的分離方法[1]。精餾塔的嚴格優化與設計對于降低精餾過程能耗，節約生產成本，減少生產過程的碳排放和環境污染具有重要意義[2-4]。

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

在上述精餾塔的優化工作中，研究者通過使進料流股和回流流股同時進入多個平衡級，建立進料位置和總塔板數優化的超結構。這種精餾塔模型的超結構較為復雜，因此很難應用到復雜精餾流程中。本文從結構優化的角度，將精餾塔優化中較為困難的進料位置和總塔板數優化轉化為空間結構的分布優化，即在精餾段和提餾段預設多于需要的平衡級數，采用結構參數代表平衡級是否存在，再通過精餾段和提餾段的實際平衡級數來表示進料位置和總塔板數。

在結構優化的研究中，信賴域優化算法更常被用來求解優化問題。不同于線搜索優化算法，信賴域優化算法采用先步長、后方向的策略，其在當前迭代點的鄰域內構造優化子問題，并逐步逼近原問題的最優解。信賴域算法具有很好的適用性和穩健性，并且可以保證優化的全局收斂性[15]，另外，信賴域算法可以方便地與流程模擬軟件進行聯用，因此，其在精餾塔優化中有許多優勢。

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

1 模型與優化算法

1.1 精餾模型

本文從結構優化的角度建立精餾模型，將總塔板數和進料位置的優化轉化為塔板空間分布的結構優化，并用0-1二元結構參數表示塔板的存在性。如圖1 所示，在精餾段和提餾段預先設置多于需要的平衡級數，每一個方框代表一個平衡級，當二元結構參數為1時，代表平衡級起到分離作用（灰色），取0 則表示平衡級不存在。這樣，通過優化0-1 二元結構參數的分布，就可以實現對總塔板數和進料位置的同時優化，將精餾過程的優化轉化為二元結構參數分布結構的優化。

圖1 精餾塔結構參數模型示意圖Fig.1 Schematic diagram of distillation column structure parameter model

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

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

1.2 移動漸近線算法

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

圖3 MMA算法Fig.3 MMA algorithm

1.3 優化方法

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

圖4 精餾塔優化方法Fig.4 Optimization method of distillation column

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

2 算例分析

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

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

在本算例中，一個單精餾塔被用于分離乙苯和苯乙烯的混合物，混合物的進料情況如圖5 所示，規定塔頂和塔釜產品的摩爾分數高于95%。因為苯乙烯在壓力較高時易于發生聚合，因此塔頂的壓力選為6 kPa。優化變量包括每一平衡級的繞流效率、回流比（RR）和再沸器的汽化分率（Vf）。本算例為每個塔節設置25 個平衡級，每個塔節有25 個繞流效率作為決策變量，因此該優化問題的決策變量數為52，約束個數為2。結合目標函數和精餾過程應滿足的約束條件，可以得到式（3）所示的約束問題，該問題的結構同樣適用于后面兩個優化算例。

圖5 乙苯-苯乙烯分離塔的結構Fig.5 Structure of ethylbenzene-styrene column

本文的優化初始值選擇回流比為5，再沸器的汽化分率為0.5，繞流效率分別從5 個不同的初始點0.1、0.3、0.5、0.7 和0.9 進行優化，不同初始點的優化結果如表1 所示。表中的N1和N2分別表示精餾段和提餾段繞流效率之和，即精餾段和提餾段的塔板數。

表1 乙苯-苯乙烯分離塔的優化結果Table 1 Optimization results of ethylbenzene-styrene column

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

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

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

2.2 直接精餾塔序列

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

圖7 直接序列精餾塔結構Fig.7 Structure of direct sequence distillation column

如表2 所示，在5 個不同初始點的優化計算中，從不同初始值開始優化，繞流效率最終均收斂為0或1，相應的平衡級數為整數，這說明在本算例中MMA 優化算法同樣可以找到局部最優解。本算例優化變量比2.1 節算例多一倍，在這種情況下，不同初始值均能收斂到局部最優解，說明MMA 優化算法的收斂性和穩定性較好。

表2 直接序列精餾塔的優化結果Table 2 Optimization results of direct sequence distillation column

2.3 部分熱耦合的精餾過程

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

圖8 部分熱耦合精餾結構Fig.8 Structure of partially thermally coupled distillation process

如表3 所示，在5 個不同初始點的優化計算中，不同初始值的優化結果的繞流效率均為0或1，這說明繞流效率模型在具有熱耦合精餾結構中依然有效。5 個不同的初始值均獲得了優化結果，這一方面說明所提優化框架穩健性較好，對于含有熱耦合結構的優化過程依然有效，另一方面說明MMA 優化算法有較好的收斂性，在含有熱耦合結構的優化中依然有效。

表3 部分熱耦合精餾過程優化結果Table 3 Optimization results of partially thermally coupled distillation process

3 結 論

（1）從結構優化的角度，將精餾塔總塔板數和進料位置的優化轉化為二元結構參數分布結構的優化，并結合繞流效率模型將二元離散變量轉化為連續變量。該模型可以降低精餾塔模型的復雜程度，實現對總塔板數和進料位置的同時優化。

（2）采用聯立方程法對精餾塔進行建模，采用虛擬瞬態連續性方程來輔助精餾穩態模擬計算，實現了精餾塔決策變量的同時優化。分別對一個單一精餾塔、直接序列精餾塔和熱耦合結構精餾塔進行優化設計。結果表明，虛擬瞬態連續性方程可有效輔助穩態模擬的計算。

（3）結合MMA 算法的特點，將Aspen Custom Modeler 和Matlab 耦合鏈接，充分利用了Aspen Custom Modeler 的模擬能力和Matlab 的優化計算能力。該優化方法可保證本文3個優化算例的正常運行，體現了該方法的穩健性。

從優化算例來看，初始值的選取會對優化結果產生影響，如何選取初始值以獲得更好的結果值得進一步探究。本文建立的優化方法對許多優化算法同樣適用，因此，不同優化算法在精餾塔優化中的表現是未來一個值得研究的方向。