基于模擬退火算法的換熱網絡雙層優化方法

彭富裕，崔國民，陳家星

（上海理工大學 新能源科學與工程研究所，上海 200093）

換熱網絡作為化工過程中重要的能量回收環節，其設計的合理性和高效性將直接關系到系統的整體性能。1969年，Masso等［1］首次嚴格定義換熱網絡綜合問題，即充分利用工藝物流的能量，使得熱回收量最大或經濟費用投入最小。以夾點技術［2］為代表的熱力學方法憑借其物理意義清晰、操作簡單的優勢得到了廣泛應用［3-5］，但此類方法分步優化的特點不能同時兼顧換熱器數、換熱面積及能量回收之間的權衡關系，只能得到接近最優的換熱網絡設計。基于數學理論提出的換熱網絡同步綜合模型由于同時涉及代表結構的整型變量和決定費用大小的連續變量優化，屬于典型的混合整數非線性規劃問題，目標函數呈現的嚴重非凸非線性使得確定性優化方法的求解過程極易陷入某個局部最小解，而無法找到理論上的全局最優解［6-7］。此外，Furman等［8］證明換熱網絡綜合問題屬NP-hard范疇，求解的可操作性嚴重受制于網絡規模的大小。

基于隨機技術的啟發式算法對目標函數的要求相對寬松，并憑借其可操作性強和計算效率高等優點開始受到越來越多的關注［9-13］。換熱網絡的優化過程既包括換熱器和公用工程面積費用以及公用工程費用的連續變量優化，同時也涉及到代表換熱器有無的0-1整型變量優化。單從整型變量角度進行分析，結構優化過程屬于大規模的組合優化問題，由于各個網絡結構之間并不存在任何聯系，因此現有的算法通常以費用或能量目標作為判斷結構優劣的標準，將進化的迭代次數或目標函數的變化情況作為最終的收斂判據。而固定結構下的連續變量優化的實質是解決如何高效跳出局部極值陷阱問題，隨著網絡規模的擴大，優化變量的數目也逐漸增加，要求算法具有較強的全局搜索能力和計算效率。Lewin等［14-15］將遺傳算法應用到換熱網絡的結構優化中，以最大能量回收為目標優化特定的網絡結構，再反饋到遺傳算法的適應函數進化網絡結構。隨后，這種采用遺傳算法優化外層結構的策略又分別與微分進化、模擬退火和粒子群等算法相結合應用于換熱網絡的最優化研究［16-19］。但隨機算法的本質卻決定了此類算法不能搜到理論上的全局最優解。

與其他隨機算法相比，模擬退火算法的一大優勢在于理論較為完善，在條件合適的前提下，其全局收斂性可以得到保證［20］。在模擬退火算法的基本思想提出之后，Kirkpatrick等［21］最早將其成功應用于組合優化問題。隨后，Dolan等［22-23］又將該算法引入到換熱網絡的設計中，用于求解同時具有離散變量和連續變量換熱網絡綜合問題，以Metropolis準則同步更新對應的換熱網絡結構和換熱量分配。在此基礎上，Athier等［24］嘗試在模擬退火機制進化外層結構的同時，采用現有非線性規劃技術優化內層連續變量，旨在提高整個算法的效率，解決工業級別的換熱網絡綜合問題。

本工作借鑒雙層優化的思想，采用換熱網絡的分級超結構模型，外層以費用目標作為判斷結構優劣的標準，將模擬退火算法應用于代表網絡結構的整型變量優化，通過特定的鄰域搜索方式及結構約束逐步進化更新網絡結構；內層費用連續變量優化同樣采用模擬退火算法進行，旨在為實現換熱網絡的全局最優化提供一種新的方法。

1 換熱網絡優化的數學模型

1.1 換熱網絡模型

以2股熱流體和3股冷流體為例，無分流的Grossmann分級超結構［3］見圖1。一般來說，換熱網絡的級數k（k=1，2，…，Nk）的取值不大于max（NH，NC）。級數確定后，不同流股間最多的匹配次數為Nk，換熱器個數最多可達到NH×NC×Nk。

1.2 目標函數及約束條件

針對由i股熱流體、j股冷流體組成的k級換熱網絡，若以經濟費用最小為目標，則優化的目標函數為：

約束條件主要包括：單股流體熱平衡、冷熱流股可行出口溫度、冷熱公用工程熱平衡及非負約束。

1）單股流體熱平衡：

2）冷熱流股可行出口溫度：

3）冷熱公用工程熱平衡：

還隱含包括以下連續變量在傳熱計算中的非負約束。

1）換熱面積非負約束：

2）換熱量非負約束：

依據模型中的假設，冷熱流股逆流布置。傳熱計算中，單個換熱器滿足以下熱平衡關系：

為了防止出現溫度交叉和換熱面積無限大的情況，必須定義相應的最小換熱溫差約束：

2 換熱網絡雙層優化策略

換熱網絡結構優化的實質是對代表換熱器有無的0-1整型變量的優化，考慮到每種整型變量的組合形式都將對應一種特定的換熱網絡結構，隨著換熱網絡規模的擴大，求解過程也變得異常復雜。采用模擬退火算法優化換熱網絡結構時，將代表結構的特定整型變量組合看作退火過程中的狀態點，通過對比不同結構的優劣決定是否接受新狀態點。對于特定結構，由于費用目標是判斷結構優劣的標準，因此要求內層的連續變量算法需具有較強的全局搜索能力，才能保證結構優劣判據的有效性。從全局最優化的角度來講，模擬退火算法允許以一定的概率接受“差”的試探點，從而可以避免結構優化和費用優化過早地陷入某個局部最小解，只要相關控制參數設置合理，可以得到問題的相對全局最優解。

2.1 內層連續變量優化

對于固定結構下的換熱網絡連續變量優化，關鍵在于如何有效地跳出局部最小解陷阱，而內層算法的計算效率也將直接影響到外層結構的進化效率，因此在保證計算精度的前提下應盡可能縮短計算時長。模擬退火算法的冷卻進度表控制著算法的全局收斂性和計算效率，其中的主要參數包括初始溫度（T0）、降溫函數、溫度迭代長度（L）及終止準則。

2.1.1 鄰域搜索方式

模擬退火算法優化換熱網絡費用連續變量，由式（1）可知，特定結構下由換熱面積及公用工程使用帶來的費用可作為優化過程中的某個狀態點。由于大量約束條件的存在，給新試探點的產生帶來一定困難，采用構建輔助函數的方法將原問題轉化為無約束問題后，也存在相懲罰因子難以確定等問題。因此，考慮將流股間匹配的換熱量（Qijk）作為優化變量，利用潛在的熱平衡約束條件，保證每次新的狀態點均在可行域中產生，具體操作如下。

1）計算各熱冷流股的換熱潛能Qpi和Qpj：

2）對于某個換熱器（zijk），給定確定換熱量增減的概率參數（β，0≤β≤1），產生［0，1］均勻分布的隨機數r和r′。

3）若r≥β，則

且

增加流股匹配換熱量。

4）若r<β，則

且

減少流股匹配換熱量。

以上鄰域搜索方法可以保證新生成的狀態點不違反能量守恒約束，均在可行域中。對于約束式（18）和（19），可以規定當出現溫度交叉時，該冷熱流股為不合理匹配，令zijk=0，作無效匹配處理。

2.1.2 初始溫度和溫度迭代長度

T0和L是影響模擬退火全局收斂性的重要因素。實際應用中，常將L設為一個常數，即在每一溫度下，均迭代相同的次數。隨后，可通過取值計算估計和統計推斷估計來確定與L相匹配的T0，具體計算過程如下：

1）給定溫度變化常量ΔT和T0的估計值、初始溫度下解的接受比率的預期值γ0（接近于1）、收斂精度ε0以及每一溫度下的L。令迭代計數k0=1。

2）記錄L步迭代后接受的狀態個數L′，該溫度下接受的比例Rk0=L′/L。

3）若─Rk0-γ0─<ε0，則停止計算，T0即為合適的初始溫度；否則，

當Rk0-1<γ0且Rk0<γ0時，k0=k0+1，T0=T0+ΔT，轉回步驟2）；

2.1.3 溫度下降函數和終止準則

溫度下降方法是影響模擬退火算法全局搜索性能的又一個重要因素。在實際應用中，往往采用較為直觀的降溫方式，即Tk+1=αTk，其中α為0到1之間的常數，每次溫度都以相同的比率下降，α越大則溫度下降越慢。理論上，當退火溫度為零時，模擬退火算法終止。但在實際應用時，往往采用一些比較直觀的終止準則，即T<ε。

2.2 外層整型變量優化

單從變量形式上看，由特定整型變量組合代表的換熱網絡結構之間并不存在任何聯系，故選取對應結構下的最優費用作為判斷結構優劣的標準。同樣，模擬退火算法用于換熱網絡的結構優化也需要確定相應的冷卻進度表：初始溫度（T0out）、降溫函數、溫度迭代步長（L0out）和終止準則。

當前結構下生成新試探結構的過程對算法的整體效率有著重要的影響，不僅要求實現整個求解域的搜索過程，還要盡可能排除明顯較“差”的換熱網絡結構以提高計算效率。因此，在生成新試探結構時，可以增加一些相應的約束條件，具體操作如下。

1）給定確定換熱器增減的概率參數βout，生成［0，1］分布的隨機變量rout。

2）若rout≥βout，從已有換熱器中隨機抽取一項（不包括公用工程），令zijk=0，刪去一個換熱器。

3）若rout<βout，從換熱網絡中隨機抽取任一不存在換熱器的位置（不包括公用工程），令zijk=1，增加一個換熱器。

新試探結構的生成過程中可設置最小總換熱器個數約束，而對單股流體也可設置最大換熱器個數約束，以防止出現換熱器總數過少或某股流體匹配過多換熱器的不合理結構。此外，對于不同的換熱網絡設計問題，通過觀察單股流體的溫位水平，可以限制“過冷”熱流體與“過熱”冷流體之間的禁忌匹配。

Lout可取固定步數，由于每個結構都對應一次費用連續變量的全局優化過程，為了提高計算效率，結構優化的初始溫度可采用估算所得，即T0out=δ(Fmax-Fmin)，δ為充分大的數。同時，溫度下降方式和終止準則沿用Tk+1out=αoutTkout和Tout<εout。

3 算例與分析

3.1 算例一

選用文獻［25］中的6×4算例，該換熱網絡由6股熱流體和4股冷流體組成，具體的物流參數詳見表1，換熱器面積費用計算公式為60A $/a，熱公用工程費用系數為100 $（kW·a），冷公用工程費用系數為15 $（kW·a）。

表1 算例一的物流參數Table 1 Problem data for case 1

選取網絡級數為4，內層連續變量優化相關參數為：β=0.95，充分利用冷熱流股換熱潛能， L=500，T0=2×106，ΔT=5×105，γ0=0.95，ε0=1×10-6，α=0.75，ε=1。外層整型變量優化相關參數為：βout=0.5，具有相同的概率增減換熱器； T0out=105， Lout=60，αout=0.9。在生成試探結構時增加以下兩個約束：

1）初始換熱器數估計值為8，規定熱流股上最大換熱器數為3、冷流股換熱器數為4（不包括公用工程）。

2）熱流股4的溫位水平較低，認為屬“過冷”熱流體，限制其不參與換熱匹配。

采用Fortran編程，計算機CPU為Intel（R）Core（TM）i5-2320，主頻3.00 GHz，計算時長為143 h。優化結果見表2，換熱網絡結構見圖2。

由表2可知，基于模擬退火算法的結構優化的實質是各個換熱網絡結構之間的相互變換，并逐步收斂于費用較低的結構的過程，可以搜索到較多的較優結構。再借助模擬退火算法在連續變量優化方面的較強全局搜索能力可以得到優于文獻結果的換熱網絡設計。此外，由圖2可知，對于熱流股4由于增加了不參與換熱匹配的約束，只能用公用工程冷卻，而單股流體的最大換熱器數約束也可以有效地將換熱網絡結構限制在合理的范圍之內。

表2 算例一的優化結果Table 2 Solutions to case 1

圖2 算例一的優化結構Fig.2 Results for case 1.

3.2 算例二

算例二選自文獻［28］，由8股熱流體和7股冷流體組成，具體物流參數見表3。與算例一的不同在于計算投資費用時考慮了固定投資費用，換熱器面積費用的計算式為8 000+500A0.75$/a，熱公用工程費用系數為80 $（kW·a），冷公用工程費用系數為10 $（kW·a）。

選取換熱網絡級數為2，內層連續變量優化相關參數為：β=0.95，L=500，T0=106，ΔT=106，γ0=0.95，ε0=1×10-6，α=0.75，ε=1。外層整型變量優化相關參數為：βout=0.5，T0out=2×106，Lout=50，αout=0.9。此外，在生成試探結構時增加以下兩個約束：

1）初始換熱器數估計值為9，規定熱流股上最大換熱器數為4、冷流股器數為3（不包括公用工程）。

2）將熱流股7與冷流股4定為固定匹配，不參與其他流體的匹配也可以各自達到目標溫度。

采用Fortran編程，計算機CPU為Intel（R）Core（TM）i5-2320，主頻3.00 GHz，計算時長為183 h。優化結果見表4，換熱網絡結構見圖3。

由表4可知，基于模擬退火算法的雙層優化方法可以優化得到優于文獻結果的換熱網絡設計。由圖3可知，雖然沒有針對冷流股7增加相應的匹配限制約束，但由于該流股屬于溫位較高的冷流體，不利于與溫位較低的熱流體進行匹配換熱，因此得到的優化結果中該流股通過公用工程加熱至目標溫度。

此外，根據歐拉廣義網絡理論對換熱網絡換熱單元個數的估計，算例一換熱單元應為11個，但實際優化得到的換熱單元分別（對應圖4中的兩個結果）為14和15個；而算例二的理論換熱單元為16個，實際換熱單元個為17和18個。設備投資費用由固定投資費用和面積費用兩部分組成。由此可知，當面積費用系數B=1且固定投資費用計算常數（CF）相對較小時，換熱網絡設計受換熱單元個數的限制也相對較小。反之，目標函數將受制于換熱網絡的換熱單元個數限制，優化過程更傾向于接受個數較少的“大面積換熱器”而不是個數眾多的“小面積換熱器”。算例一的固定投資費用為0、B=1，因此受換熱單元個數的限制較小，得到的優化結果換熱單元個數均遠大于理論值。而算例二由于存在較大的固定投資費用且B=0.75，因此最優的換熱網絡結構更傾向于換熱器更少的結構。

表3 算例二的物流參數Table 3 Problem data for case 2

表4 算例二的優化結果Table 4 Solutions to case 2

圖3 算例二的優化結構Fig.3 Results for case 2.

3.3 討論

效率和精度是評價優化方法的兩個重要指標。雙層優化的思想將同時具有整型變量和連續變量的混合整數非線性規劃問題分為外層整型變量優化和內層連續變量優化兩個過程。相對于整型變量的全局最優化，針對連續變量的全局最優化技術相對成熟，從計算精度角度來看，無論是采用等效松弛處理非凸項的確定性方法還是眾多隨機性方法，都可以有效實現連續變量的全局最優。整型變量的全局最優化一直以來都是優化理論中的難點問題，代表換熱網絡結構的整型變量優化，其實質是復雜的組合優化問題，對于由NH股熱流體和NC股冷流體組成的Nk級換熱網絡，理論上存在2NH×NC×Nk種可能的結構組合，而隨著換熱網絡規模的擴大和優化變量的增加，易導致組合爆炸問題的發生，從而大幅降低問題的可操作性，進而影響優化的效率。

因此，為了同時兼顧雙層優化方法的效率和精度，要保證外層算法在有限的時間內盡可能試探更多的結構組。本工作從傳熱角度添加額外結構約束的方法旨在縮小搜索域，減少對部分較差結構的計算時間。但由于模擬退火算法中新試探點是在當前解的鄰域范圍內通過隨機擾動產生的，因此這種方法也可能會限制結構產生過程中結構與結構之間的過渡過程。今后的研究可致力于更有效的結構約束的確定。為了提升整個算法的計算效率，內層連續變量優化要在保證全局收斂性的前提下，盡可能地縮短單個結構的優化計算時間，因此可以考慮采用更高效的連續變量全局最優化方法。

4 結論

1）針對同時存在整型變量和連續變量的換熱網絡優化問題，提出一種基于模擬退火算法的雙層優化方法，分別處理兩種變量形式的優化過程。

2）外層算法優化選取費用目標作為判斷結構優劣的標準，通過隨機擾動的方式不斷產生新的結構，并借助模擬退火機制不斷進化換熱網絡結構；內層算法通過模擬退火算法對固定結構下的連續變量進行優化，再將得到的最優解反饋到外層算法，作為結構優化的判斷標準。

3）通過兩個具體算例證明，模擬退火算法可有效處理整型變量形式的結構優化，找到較優結構，而雙層優化方法可以有效處理具有整型變量和連續變量兩種變量形式的換熱優化問題，取得較好的換熱網絡設計。

符 號 說 明

A 換熱器面積，m2

B 投資費用計算常數

C 投資和運行費用計算常數

CF固定投資費用計算常數

Fmax，Fmin目標函數的最大值和最小值

Fcp熱容流率，kW/℃

k0初始溫度計算迭代計數

L 溫度迭代步長

Q 換熱量，kW

Qp換熱潛能，kW

Rk0解的接受比例

r ［0，1］內均勻分布的隨機數

T 溫度

T0初始溫度

ΔT 溫度變化常量

ΔTmin最小傳熱溫差

ΔTLM對數平均溫差

U 換熱系數，kW/（m2·℃）

z 換熱器有無的0-1變量

α 溫度下降系數

β 概率參數

γ0初始溫度下解的接受比率的預期值

ε 終止溫度

ε0收斂精度

上角標

IN 流股入口

in 單個換熱器入口

OUT 流股出口

out 單個換熱器出口

k 迭代次數

下角標

CU 冷公用工程

HU 熱公用工程

i 熱流股，i=1，2，…，NH

j 冷流股，j=1，2，…，NC

k 換熱網絡級數，k=1，2，…，Nk

out 外層算法

［1］Masso A H，Rudd D F.Synthesis of System Designs：Ⅱ.Heuristic Structuring［J］.AIChE J，1969，15（1）：10-17.

［2］Linnhoff B，Flower J R.Synthesis of Heat Exchanger Networks：Systematic Generation of Energy Optimal Networks［J］.AIChE J，1978，24（4）：633-654.

［3］Rivedi K K，Oeill B K，Roach J R，et al.A New Dual-Temperature Design Method for the Synthesis of Heat Exchanger Networks［J］.Comput Chem Eng，1989，16（6）：667-685.

［4］Fraser D M.The Use of Minimum Flux Instead of Minimum Approach Temperature as a Design Specification for Heat Exchanger Networks［J］.Chem Eng Sci，1989，44（5）：1121-1127.

［5］Trivedi K K，Oeill B K，Roach J R，et al.Systematic Energy Relaxation in MER Heat Exchanger Networks［J］.Comput Chem Eng，1990，14（6）：601-611.

［6］Yee T F，Grossmann I E.Simultaneous Optimization Models for Heat Integration：Ⅱ.Heat Exchanger Network Synthesis［J］.Comput Chem Eng，1990，14（10）：1165-1184.

［7］胡向柏，崔國民，涂惟民.復雜換熱網絡MINLP中的非線性特性分析［J］.工程熱物理學報，2012，33（2）：285-287.

［8］Furman K C，Sahinidis N V.Computational Complexity of Heat Exchanger Network Synthesis［J］.Comput Chem Eng，2001，25（9/10）：1371-1390.

［9］Chakraborty S，Ghosh P.Heat Exchanger Network Sythesis：The Possibility of Randomization［J］.Chem Eng J，1999，72（3）：209-216.

［10］Pariyani A，Gupta A，Ghosh P.Design of Heat Exchanger Networks Using Randomized Algorithm［J］.Comput Chem Eng，2006，30（6/7）：1046-1053.

［11］Gupta A，Ghosh P.A Randomized Algorithm for the Ef fi cient Synthesis of Heat Exchanger Networks［J］.Comput Chem Eng，2010，34（10）：1632-1639.

［12］Lin B，Miller D C.Solving Heat Exchanger Network Synthesis Problems with Tabu Search［J］.Comput Chem Eng，2004，28（8）：1451-1464.

［13］Silva A P，Ravagnani M A S S，Biscaia E C，et al.Optimal Heat Exchanger Network Synthesis Using Particle Swarm Optimization［J］.Optim Eng，2010，11（3）：459-470.

［14］Lewin D R，Wang H，Shalev O.A Generalized Method for HEN Synthesis Using Stochastic Optimization：Ⅰ.General Framework and MER Optimal Synthesis［J］.Comput Chem Eng，1998，22（10）：1503-1513.

［15］Lewin D R.A Generalized Method for HEN Synthesis Using Stochastic Optimization：Ⅱ.The Synthesis of Cost-Optimal Networks［J］.Comput Chem Eng，1998，22（10）：1387-1403.

［16］Yerramsetty K M，Murty C V S.Synthesis of Cost-Optimal Heat Exchanger Networks Using Differential Evolution［J］.Comput Chem Eng，2008，32（8）：1861-1876.

［17］Khorasany R M，Fesanghary M.A Novel Approach for Synthesis of Cost-Optimal Heat Exchanger Networks［J］.Comput Chem Eng，2009，33（8）：1363-1370.

［18］Luo Xing，Wen Qingyun，Fieg Georg.A Hybrid Genetic Algorithm for Synthesis of Heat Exchanger Networks［J］.Comput Chem Eng，2009，33（6）：1169-1181.

［19］Huo Zhaoyi，Zhao Liang，Yin Hongchao，et al.Simultaneous Synthesis of Structural-Constrained Heat Exchanger Networks with and Without Stream Splits［J］.Can J Chem Eng，2013，91（5）：830-842.

［20］Dekkers A，Aarts E.Global Optimization and Simulated Annealing［J］.Mathematical Programming，1991，50（1/3）：367-393.

［21］Kirkpatrick S，Gelatt C D，Vecchi M P.Optimization by Simulated Annealing［J］.Science，1983，220：671-680.

［22］Dolan W B，Cummings P T，LeVan M D.Process Optimization Via Simulated Annealing：Application to Network Design［J］.AIChE J，1989，35（5）：725-736.

［23］Dolan W B，Cummings P T，LeVan M D.Algorithmic Ef fi ciency of Simulated Annealing for Heat Exchanger Network Design［J］.Comput Chem Eng，1990，14（10）：1039-1050.

［24］Athier G，Floquet P，Pibouleau L，et al.Synthesis of Heat-Exchanger Network by Simulated Annealing and NLP Procedures［J］.AIChE J，1997，43（11）：3007-3020.

［25］Ahmad S.Heat Exchanger Networks：Cost Trade-Offs in Energy and Capital［D］.Manchester：Manchester Institute of Science and Technology，1985.

［26］Ravagnani M A S S，Silva A P，Arroyo P A，et al.Heat Exchanger Network Synthesis and Optimization Using Genetic Algorithm［J］.Appl Therm Eng，2005，25（7）：1003-1017.

［27］萬義群，崔國民，彭富裕，等.基于溫差原則的換熱網絡分步綜合策略［J］.石油化工，2013，32（9）：996-1003.

［28］Bjork K M，Pettersson F.Optimization of Large-Scale Heat Exchanger Network Synthesis Problems［C］//Hamza M H，ed.Palm Springs：Proceedings of the IASTED International Conference on Modeling and Simulation，2003：313-318.

［29］胡向柏，崔國民，涂惟民，等.網絡填充函數法的全局優化［J］.化學工程，2011，39（1）：28-31.

［30］張勤，崔國民，關欣.基于蒙特卡羅遺傳算法的換熱網絡優化問題［J］.石油機械，2007，35（5）：19-22.

［31］胡向柏，崔國民，許海珠，等.換熱器優化順序對換熱網絡全局優化的影響［J］.化工進展，2012，31（5）：987-1003.