基于分子熱力學(xué)性質(zhì)和密度峰聚類的脫硫汽油集總

李懷旭，孫曉巖，陶少輝，夏力，項曙光

（青島科技大學(xué)化工學(xué)院，過程系統(tǒng)工程研究所，山東 青島 266042）

引 言

隨著石油加工領(lǐng)域進(jìn)入分子管理時代，業(yè)內(nèi)逐漸形成一套完整的油品分析模擬體系，即油品分子表征（重構(gòu)）-組分集總（劃分）-解集總（還原）體系。油品分子表征通過現(xiàn)代化分析儀器進(jìn)行[1-3]，在此基礎(chǔ)上依靠計算機(jī)輔助方式構(gòu)建分子體系的模型化方法來進(jìn)行分子重構(gòu)[4-7]。油品真實分子組成在帶來更豐富信息的同時，也極大增加了煉油過程模擬和優(yōu)化模型的規(guī)模，導(dǎo)致計算量激增，因而如何在盡量保留油品真實分子組成信息的同時，有效減少組分?jǐn)?shù)以提高煉油過程的模擬與優(yōu)化效率，是石油分子管理體系亟待解決的關(guān)鍵問題之一。而集總過程則可以根據(jù)一定規(guī)則將油品分子劃分為少數(shù)幾個虛擬組分，是減小過程模擬計算規(guī)模的有效方式[8-9]。如在油藏模擬領(lǐng)域，Montel 等[10]首次提出通過迭代的動態(tài)聚類（dynamical clustering,DynC）算法將大量組分通過性質(zhì)的相似性進(jìn)行集總，將150 個組分用少于10 種虛擬組分替代，并提出了虛擬組分物性的計算方法?；谏倭刻摂M組分進(jìn)行模擬后，為獲得模擬結(jié)果中的產(chǎn)品真實組成，Leibovici 等[11]通過閃蒸計算實現(xiàn)解集總（delumping）來還原油品組分，封松等[12-13]在其構(gòu)建的分子級組成模型中利用可逆集總方案的映射規(guī)則倒推產(chǎn)品分子組成。

Montel 等[10]對真實組分進(jìn)行集總時，采用了需要事先指定類別數(shù)和初始聚類中心的迭代動態(tài)聚類算法，但針對特定體系，確定虛擬組分?jǐn)?shù)目需要進(jìn)行大量的分析和比較，且聚類結(jié)果對于初始中心的選擇極為敏感[14-16]。這些弊端使得動態(tài)聚類在實際應(yīng)用時難以產(chǎn)生穩(wěn)定的結(jié)果，即使Montel 等提出通過人工選取類別數(shù)量最多的組分為中心可以使算法變得更穩(wěn)定，但這是基于經(jīng)驗的，依賴于使用者對數(shù)據(jù)集的了解程度。對于過程模擬計算而言，往往希望得到一套穩(wěn)健高效的集總計算模型，以適用于不同類型的數(shù)據(jù)。

基于此，本工作對油品分子的虛擬組分劃分?jǐn)?shù)目（即集總數(shù)目）對模擬計算的影響進(jìn)行討論分析，將機(jī)器學(xué)習(xí)領(lǐng)域中一種無須事先指定類別數(shù)和聚類中心的密度峰聚類（density peak clustering,DPC）[17]算法引入油品分子的集總過程，以真實分子的熱力學(xué)性質(zhì)為聚類變量，從而實現(xiàn)對表征汽油分子的集總，同時以相平衡與精餾計算對方法的有效性進(jìn)行驗證。

1 模型與方法

1.1 DPC算法原理

DPC 算法基于一個核心假設(shè)，即認(rèn)為給定數(shù)據(jù)集的中心一定是該數(shù)據(jù)集中密度較大且互相遠(yuǎn)離（相對距離較大）的一些數(shù)據(jù)點。因而該算法的計算過程取決于數(shù)據(jù)集的距離矩陣dn×n（n為樣本個數(shù)）和每個點的密度ρi。距離是聚類算法中常用的相似性度量標(biāo)準(zhǔn)，采用何種距離計算公式取決于不同種類的算法和目標(biāo)數(shù)據(jù)集的分布情況。這里采用Montel等[10]提出的距離公式來度量數(shù)據(jù)點i與j間的差異，如式（1）所示。

式中，wk為不同性質(zhì)（維度）的權(quán)重值；pik與pjk為數(shù)據(jù)點的性質(zhì)值；N為選取性質(zhì)的個數(shù)（即數(shù)據(jù)集的維數(shù)）。采用式（1）度量性質(zhì)之間的相似性旨在擴(kuò)大性質(zhì)間的微小差距。對于數(shù)據(jù)點的密度，本文采用高斯核計算數(shù)據(jù)點的局部密度[18-19]：

式中，dc為截斷距離，其物理意義是用于判定數(shù)據(jù)點的疏密程度，定義為將距離矩陣的元素升序排列，其1%～5%處對應(yīng)的距離。需要指出，dc為該算法唯一需要人工指定的參數(shù)，Rodriguez 等[17]的研究結(jié)果表明，其在一定范圍內(nèi)調(diào)節(jié)有助于更好地發(fā)現(xiàn)中心點。

相對距離δi定義為數(shù)據(jù)點i與所有密度大于點i的數(shù)據(jù)點之間的最小距離，對于密度最高的數(shù)據(jù)點，相對距離表達(dá)式為

根據(jù)各數(shù)據(jù)點的ρ和δ可做出二維平面決策圖，在圖中明顯離群的數(shù)據(jù)點即為類中心點，后文實例應(yīng)用時將會詳細(xì)展示如何尋找中心點。得到中心點后，將其余點分配給比其密度高的最近鄰點所在的類別即可。DPC算法步驟如下所示。

（1） 對數(shù)據(jù)集進(jìn)行標(biāo)準(zhǔn)化處理，將所有性質(zhì)變量變換到[0,1]區(qū)間；

（2） 通過式（1）計算距離矩陣dn×n：

（3） 通過式（2）計算每個數(shù)據(jù)點的密度；

（4） 通過式（3）、式（4）計算所有數(shù)據(jù)點的δ；

（5） 做出決策圖并選取類中心；

（6） 按照分配機(jī)制分配其他數(shù)據(jù)點。

由上述步驟可知，DPC 算法應(yīng)用到油品分子的集總劃分時的優(yōu)點為：（1）無須預(yù)定中心點，這使得在確定虛擬組分?jǐn)?shù)目時有足夠的參考；（2）算法所得到的類中心是隸屬于數(shù)據(jù)集的真實樣本，且將每一類看作是數(shù)據(jù)空間中被低密度區(qū)域分割開的稠密區(qū)域，故能檢測到任意形狀的數(shù)據(jù)集，對于非球形數(shù)據(jù)集有很好的適用性，一些研究已經(jīng)對此加以證明[20-25]；（3）DPC 算法流程中無迭代過程，類中心的發(fā)現(xiàn)和其余點的分配都是一次性完成的。這意味著在預(yù)先計算好距離矩陣后，DPC 算法效率不受數(shù)據(jù)集維度影響，而動態(tài)聚類法雖然在處理大量數(shù)據(jù)時有良好的伸縮性，但需要反復(fù)地迭代計算距離并重置中心點，且迭代結(jié)果受初始中心的影響很大，每次的聚類結(jié)果也都存在差異[26-28]。

1.2 集總性質(zhì)選擇

Montel 等[10]提出以Peng-Robinson（PR）狀態(tài)方程參數(shù)a、b、m以及分子量四個參數(shù)作為集總性質(zhì)，其中參數(shù)a、b、m分別通過式(5)～式(7)求得。

可見，參數(shù)a和b均由臨界壓力Pc和臨界溫度Tc求出，m由偏心因子ω求解，四個參數(shù)除了能很好地反映純物質(zhì)的臨界特征和極性特征外，還能體現(xiàn)各個組分在熱力學(xué)過程中的性質(zhì)變化，且與進(jìn)行模擬計算時選取的熱力學(xué)計算方法保持一致。此外，Montel 等[10]建議a、b、m以及分子量（MW）在計算相似度時的權(quán)重值分別為1.0、1.0、0.1、0.5，這里也采用此權(quán)重值。

1.3 虛擬組分物性估算

對于集總得到的虛擬組分的物性計算，Leibovici 等[29]提出了集總系統(tǒng)特性估計方程對集總組分信息損失進(jìn)行最小化處理，該方法準(zhǔn)確且與使用的熱力學(xué)模型完全一致。但在集總組分過多的情況下，該方法計算過程過于復(fù)雜，且本文側(cè)重于研究集總劃分方法而非熱力學(xué)計算，故在虛擬組分物性計算方面，本文采用一組簡單的石油餾分物性估算模型集合來進(jìn)行計算，具體如下所述。

臨界溫度和臨界壓力采用Riazi-Daubert 模型計算，模型關(guān)聯(lián)式如式（8）、式（9）所示[30]。

式中，Tb為正常沸點；Sg為相對密度。

偏心因子采用Lee-Kesler 模型計算，關(guān)聯(lián)式如式（10）所示[31]。

式中，xi表示第i個組分；c為組分總數(shù)。相對密度（Sg）可以用分子量和沸點通過關(guān)聯(lián)式估算，故這里不再收集每個組分的密度，其他諸如液相摩爾體積、水溶性、黏度、汽化潛熱等性質(zhì)的計算模型見文獻(xiàn)[32]。

2 實例應(yīng)用與討論

2.1 數(shù)據(jù)來源

油品數(shù)據(jù)來源于某催化汽油吸附脫硫工藝體系的287 種真實組分，其中包括C2～C12烯烴、C3～C12正構(gòu)烷烴、C4～C10異構(gòu)烷烴、C5～C10環(huán)烷烴、C6～C12芳烴以及二十余種含硫烴類化合物。

為方便比較，以a、b、m和MW 四個參數(shù)為數(shù)據(jù)特征。圖1 顯示了287 個數(shù)據(jù)點的分布情況，每一個子圖顯示了兩個性質(zhì)間的相互關(guān)系。從對角線的直方圖和其余散點圖可以大致判斷出由油品分子性質(zhì)組成的數(shù)據(jù)集有三個主要特征：(1)a、b、m和MW 歸一化數(shù)值都大量集中在0.25～0.75 之間；（2）組分兩兩性質(zhì)之間的線性相關(guān)性較差（擬合優(yōu)度在0.79～0.95 間）；(3) 數(shù)據(jù)集在二維平面上呈現(xiàn)出不規(guī)則的條狀或橢圓狀形態(tài)，大部分?jǐn)?shù)據(jù)點處于高度密集且重疊的狀態(tài)，不同類別之間沒有明顯的界限。顯然這對于一般的聚類方法而言很不友好。

圖1 標(biāo)準(zhǔn)化后的數(shù)據(jù)分布Fig.1 Standardized data distribution

2.2 集總數(shù)目對計算結(jié)果的影響

對于動態(tài)聚類而言，集總數(shù)目的選擇會影響最終的集總效果，因此下文將通過指定壓力（0.1 MPa）和汽化分率以計算物流溫度的閃蒸過程，討論不同集總數(shù)目對計算結(jié)果的影響。

圖2 展示了使用動態(tài)聚類的集總方法將287 種真實組分分別集總為4、6、8、16、21種虛擬組分的閃蒸計算結(jié)果。顯然，計算結(jié)果的準(zhǔn)確性與集總數(shù)目密切相關(guān)。當(dāng)集總數(shù)目為4 個時（4 components, 4 Comp），虛擬組分的閃蒸曲線的斜率小于真實組分，因為過少的集總數(shù)目會使得每個虛擬組分中包含許多性質(zhì)差異較大的組分，從而使其熱力學(xué)狀態(tài)變化難以描述實際情況，同樣6 Comp 時亦如此。再者，計算誤差也會因系統(tǒng)中虛擬組分?jǐn)?shù)量太少而擴(kuò)大。隨著集總數(shù)目的增加，真實組分得到更細(xì)致的分配，其閃蒸溫度曲線的斜率更接近實際情況，但當(dāng)集總數(shù)目太多時（21 Comp），許多本該屬于同一類別的組分被強(qiáng)行分配到不同類別，特別是對于密集分布的數(shù)據(jù)集，這會導(dǎo)致集總類別之間出現(xiàn)鏈?zhǔn)叫?yīng)，即相鄰的每個虛擬組分的性質(zhì)都會產(chǎn)生較大的誤差。

圖2 不同集總數(shù)目在0.1 MPa下的閃蒸計算結(jié)果Fig.2 Flash calculation results at 0.1 MPa for different lumping

考慮多個集總數(shù)目下，不同汽化分率時的閃蒸溫度的平均絕對誤差（mean absolute deviation,MAE），得到如圖3 所示的變化趨勢（對0 個集總數(shù)目下的MAE做了縮小處理）。這說明集總結(jié)果的準(zhǔn)確程度并不是隨著集總數(shù)目的增加持續(xù)上升，這中間存在或大或小的波動。理論上，集總數(shù)目越多，計算規(guī)模便越大，復(fù)雜度自然也越大。本例中，集總數(shù)目在10～16 之間時，MAE 較小，集總數(shù)目持續(xù)上升至21時，MAE出現(xiàn)激增，而后又出現(xiàn)驟降。

圖3 閃蒸溫度的平均絕對誤差和計算復(fù)雜度隨集總數(shù)目增加的變化Fig.3 Variation of MAE and computational complexity of the flash temperature with the increase of the number of lumps

一般地，以誤差的大小衡量集總結(jié)果的準(zhǔn)確性，得到圖4（a）所示的變化趨勢。Opt（optimal）代表最優(yōu)集總數(shù)目，此時較少的虛擬組分?jǐn)?shù)量卻有較高的準(zhǔn)確性；Lim（limit）代表極限集總數(shù)目，此時的集總數(shù)目等于組分個數(shù)，相當(dāng)于每一個組分作為一個集總。小于Opt的部分代表許多性質(zhì)差異大的組合更多地被集總在一起。Opt 到Lim 之間代表最優(yōu)集總數(shù)目下的虛擬組分被拆開，其子組分先是被互相劃分到與其性質(zhì)差異較大的類別中，后又得到了更精細(xì)的組合，此時的復(fù)雜度也更大。圖4（b）描述了一種簡化情況來表示Opt-Lim 間虛擬組分的劃分情況，P1,P2,P3 為Opt 中的某三個虛擬組分，實線圈起的部分代表原本正確的劃分被拆開，虛線圈起的部分代表被拆開的子組分重新組合為新的虛擬組分（假設(shè)新的虛擬組分比之前劃分的更加精細(xì)）。

圖4 計算精度與組分分配情況隨集總數(shù)目變化的示意圖Fig.4 Schematic diagram of the variation of calculation accuracy and component distribution with the number of lumps M, S—增加的集總數(shù)目; v, r, z—真實子組分

此發(fā)現(xiàn)具有一定意義，說明在組分集總過程中，確實存在一種兼顧效率和準(zhǔn)確性的可能，而DPC 集總方法的作用就在于自適應(yīng)地發(fā)現(xiàn)圖中所示的Opt。需要指出，當(dāng)物性估算模型依賴于摩爾分?jǐn)?shù)權(quán)重時，這一規(guī)律并不明顯。因為對于含量很少的組分，即使被錯誤分配，對虛擬組分的性質(zhì)影響也不是很大。但對于含量高的組分（尤其是重組分），差距將會很明顯。

2.3 DPC集總過程

首先根據(jù)式（1）、式（2）計算得到距離矩陣d287×287和局部密度ρi，其中dc取1%計算（經(jīng)試驗，本例dc在1%～5%內(nèi)變化對最終聚類結(jié)果影響很?。?，再根據(jù)式（3）、式（4）得到δi，得到如圖5（a）所示的平面決策圖。圖中彩色的數(shù)據(jù)點為明顯離群的點，根據(jù)算法假設(shè)，這些數(shù)據(jù)點的共同特征是同時具有較大的δi和ρi，這意味著它們彼此相互遠(yuǎn)離且擁有各自的稠密區(qū)域，故將此視為數(shù)據(jù)集的10 個中心點。為了更清晰地展示類中心點與其他數(shù)據(jù)點區(qū)別，將每個數(shù)據(jù)點的ρ2和δ相乘，記作γ并降序排列，采用Lv 等[33]提出的一種通用閾值計算方法計算γ閾值（詳細(xì)算法過程參見該文獻(xiàn)），如圖5（b）所示，從決策圖中獲取的10 個中心點在圖5（b）中依然得到了突出顯示。就此將油品中含有的287個烴類分子劃分為10類，即10個虛擬組分。

圖5 集總數(shù)目及中心組分的決策圖Fig.5 Decision diagram of the number of lumps and center components

此外，決策圖中的黑色菱形標(biāo)記也是數(shù)據(jù)集中明顯離群的數(shù)據(jù)點，但與彩色的中心點不同，菱形點只有較大的δ，其ρ卻異常?。ń咏?），這說明此類數(shù)據(jù)點彼此孤立且相互遠(yuǎn)離，在整個數(shù)據(jù)集中不屬于任何一類，本文將其定義為“離群組分”。對于其他ρ極小的點（緊靠縱坐標(biāo)軸），因為其δ并不突出，故不能作為離群組分，這類點通常處于“團(tuán)簇”的邊緣。通過查找，3 個離群組分分別對應(yīng)正十二烷、碳二苯并噻吩、1,4-二異丙苯（按δ值降序），其部分性質(zhì)已經(jīng)列于表1，與其他組分比較發(fā)現(xiàn)，它們具有個別較為突出的性質(zhì)，因此將其視作“離群”是合理的。對于離群組分，本文的處理方式與算法要求不同（Rodriguez 等[17]認(rèn)為每個離群點應(yīng)單獨為一類），本文建議兩種處理方法：若離群組分?jǐn)?shù)量太多，則將其歸為一類，作為一個虛擬組分；若其數(shù)量較少、含量低，則將其包含在較近的類別中。本例中僅含2個離群組分，且為重組分，故采取第二種處理方法。

表1 離群組分的部分性質(zhì)Table 1 Partial properties of outlier components

需要指出，本例中含有硫化氫氣體和二氧化硫氣體，其理論上也為離群組分，但算法并未偵測到兩者，在此有必要做出解釋：DPC 集總過程使用了式（1）來度量組分間的相似性，式（1）雖能擴(kuò)大性質(zhì)間的差異，但有一個潛在問題，即對于較小的數(shù)值并不靈敏。如圖6 所示，點Ⅰ、Ⅲ代表兩個離群組分，點Ⅱ代表某正常組分，且假設(shè)A、B、C為三個組分的性質(zhì)。其中點Ⅰ的性質(zhì)比點Ⅱ低50%（反向差距），點Ⅲ的性質(zhì)比點Ⅱ高50%（正向差距），若采用歐氏距離計算空間中的真實距離，則兩者與點Ⅱ的相似度是一致的，兩者都應(yīng)該計算出較小的ρ，進(jìn)而得到較大的δ，最后形成孤立點。但當(dāng)采用式（1）時，反向差距大的組分得不到較大的相對距離，即使它看起來與正向差距大的組分并無區(qū)別，這便導(dǎo)致了所有類似的組分不能在決策圖中突出顯示。本例中，H2S 和SO2具有很小的分子量，很低的臨界溫度和很大的臨界壓力，故其a、b、m、MW 均很小，算法沒有對其進(jìn)行挖掘。實際上，諸如SO2等不凝氣或一些雜質(zhì)氣體往往單獨從油品分子中劃分出來，所以在檢測離群組分時可以不考慮這一類反向差距很大的組分；相反，重組分的偵測更加重要，同時也更難偵測，而式（1）對于正向差距大的組分敏感性更強(qiáng)，故這種處理方法是合理的。

圖7（a）為最終集總結(jié)果，圖中的不同著色代表不同的類別，即反映不同真實組分所屬的虛擬組分，對角線的直方圖匯總了每種虛擬組分的性質(zhì)分布情況。通過觀察圖1 發(fā)現(xiàn)，參數(shù)a、MW 相關(guān)性較其他更高，為更直觀地顯示集總結(jié)果，暫時忽略性質(zhì)a，僅對b、m、MW 進(jìn)行三維可視化，如圖7（b）所示。經(jīng)集總得到的虛擬組分的部分物性由式（8）～式（14）計算，與其中心組分的性質(zhì)一并列于表2。

表2 中心組分與對應(yīng)的虛擬組分的部分性質(zhì)Table 2 Partial properties of center components and corresponding pseudo-components

圖8展示了每個虛擬組分中子組分的分子類型和碳原子數(shù)?？梢?，PC4 和PC5 分別主要包含C8異構(gòu)烷、C10異構(gòu)烷，兩者數(shù)量較少，即所謂的“小簇”。關(guān)注數(shù)據(jù)集中的“小簇”極為重要，這在一定程度上能反映算法對數(shù)據(jù)進(jìn)行挖掘的徹底性。從圖7（b）中又可以發(fā)現(xiàn)兩個“小簇”都分布在PC6和PC10中，而PC6 包含較多C8～C9的異構(gòu)烷及環(huán)烷，PC10 則主要為C11以上的芳烴。

圖7 287個組分的集總結(jié)果分布及結(jié)果可視化Fig.7 Lumped results of 287 components and visualization of results

表面上，不管從數(shù)據(jù)集分布形態(tài)還是具體組分來看，似乎PC4 應(yīng)隸屬于PC6。但實際發(fā)現(xiàn)，PC4 中10 個組分的平均偏心因子為0.38，而PC6 中39 個組分的平均偏心因子僅為0.30，且有近半數(shù)的組分的偏心因子小于0.30，這說明PC6 中組分的極性與PC4 有明顯差距，將此團(tuán)數(shù)據(jù)集分成兩類是符合實際的。對于PC5 而言，結(jié)合圖8（a）和圖8（b）可以很直觀地發(fā)現(xiàn)其與PC10 中的組分性質(zhì)差距更大，PC10 以C11以上的芳烴為主，其極性較大，而PC5 以C10以下的異構(gòu)烷烴為主，其極性要小于芳烴。再者，通過比較發(fā)現(xiàn)，PC5 中組分的平均臨界壓力為2.17 MPa，平均臨界溫度為337℃，而PC10 為2.52 MPa、415℃，兩者存在明顯差距。而DPC 集總方法恰好可以有效地對PC4 和PC5 進(jìn)行挖掘，這得益于該算法擅長處理具有重疊性的非球形數(shù)據(jù)集。

圖8 虛擬組分中子組分的烴類分子構(gòu)成和碳原子數(shù)分布Fig.8 Hydrocarbon composition and carbon number distribution of the subcomponents of pseudo-component

2.4 相平衡計算驗證

本文重在考察虛擬組分在汽液相平衡計算中的表現(xiàn)。Leibovici 等[11]在其體系已證實了基于動態(tài)聚類的集總方法在閃蒸計算中的可靠性。需要指出，基于動態(tài)聚類的油品集總方案需要預(yù)先指定集總數(shù)目，為便于比較，將基于DPC 算法集總得到的10 個虛擬組分作為其預(yù)先設(shè)定的集總數(shù)目。下面展示了DPC 集總方法和動態(tài)聚類方法與真實組分的對比。

很明顯，DPC 方法取得了預(yù)期效果。圖9 中，其閃蒸溫度曲線與真實組分的貼合程度較高；而采用動態(tài)聚類的集總方法即使預(yù)先給定了一個近乎合理的集總數(shù)目，其效果與DPC 方法相比也有一些差距，尤其是在接近露點溫度時誤差較大，原因是動態(tài)聚類無法捕捉不規(guī)則或具有層疊性的數(shù)據(jù)，故部分?jǐn)?shù)據(jù)點被劃分到了與自身性質(zhì)差異較大的類別中，進(jìn)而在宏觀上引起了計算結(jié)果的差異。從圖9得出，DPC 集總方法所得到的閃蒸溫度的MAE 較動態(tài)聚類方法而言下降36.0%，計算精度有明顯提升。

圖9 不同集總方法在0.1 MPa下的閃蒸溫度曲線對比Fig.9 Comparison of flash temperature curves at 0.1 MPa for different lumped methods

圖10為不同集總方法在0.1 MPa下的閃蒸熱負(fù)荷的絕對誤差隨汽化分率的變化。與真實情況相比，DPC 方法得到的各點誤差均在較小范圍之內(nèi)，而基于動態(tài)聚類的方法只有在汽化分率0.7～0.8 之間才能達(dá)到與其持平的效果。從圖10得出，DPC 集總方法所得到的閃蒸熱負(fù)荷的MAE 較動態(tài)聚類方法而言下降64.2%，計算精度有明顯提升。熱負(fù)荷是對虛擬組分的熱容、焓值等性質(zhì)的綜合考察，故此結(jié)果在一定程度上說明了DPC 集總方法的可靠性。

圖10 不同集總方法下閃蒸熱負(fù)荷的絕對誤差（與真實組分下的熱負(fù)荷相比）比較Fig.10 Comparison of the absolute error of the flash heat duty(compared to the heat load of the real component) of different lumped methods

2.5 中心組分和虛擬組分的比較

前文已述，通過DPC 算法得到聚類中心隸屬于真實組分，其在一定程度上能代表子組分的性質(zhì)，故有必要討論中心組分是否可以表征油品分子的真實性質(zhì)。圖11、圖12展示了中心組分和虛擬組分分別表征油品時的結(jié)果對比。由圖可見，在汽化分率小于0.1 時，中心組分代替真實組分的效果尚好，但隨著汽化分率的增大，其閃蒸溫度總是低于真實情況，這說明中心組分中有過多的輕組分（相比于虛擬組分而言較輕）。從圖12 的PT 包絡(luò)圖來看，臨界條件的低估又說明中心組分的極性和沸點難以反映真實情況。表3中列出了每個虛擬組分與中心組分的沸點差距，可見PC6 的中心組分的沸點比虛擬組分低28.33 K，PC8 的中心組分的沸點比虛擬組分低32.60 K。

圖11 中心組分與虛擬組分在不同汽化分率下的閃蒸溫度對比（0.1 MPa）Fig.11 Comparison of the flash temperature of the center component and the pseudo-component at different vaporization fractions at 0.1 MPa

圖12 中心組分與虛擬組分的PT相包絡(luò)圖對比（1 bar=0.1 MPa）Fig.12 Comparison of the PT phase envelopes of center component and pseudo-component

表3 各虛擬組分與其中心組分的沸點差距Table 3 Boiling point gap between each pseudocomponent and its center component

雖然中心組分目前不能完全反映虛擬組分的性質(zhì)，但分析其中的原因尚有價值。實際上，通過算法得到的虛擬組分，可看作是空間中一塊稠密程度從中心向周圍不規(guī)則發(fā)散的區(qū)域。如圖13（a）所示，以b和m為特征性質(zhì)，將PC8 中子組分的局部密度ρ映射到氣泡圖中，圖中面積最大的點為中心組分，可見，ρ較大的組分全部靠近中心組分，這些組分的性質(zhì)與中心組分更相似；而ρ較小的組分距離中心組分較遠(yuǎn)，其與中心組分的相似度也就較差（相比于本類別中的ρ較大的組分而言）。更直觀地，圖13（b）的密度分布熱力圖展示了中心組分對其子組分的描述能力大小。熱力圖的核心部分（黑色區(qū)域）為中心組分，它對其周圍的密度接近自己的點能很好反映；但對稍微遠(yuǎn)離的點（綠色及藍(lán)色區(qū)域），其代表性是有限的。對于包含子組分較多的虛擬組分，由于當(dāng)前數(shù)據(jù)分布的不規(guī)則性，其中心組分很難描述到每一個組分的性質(zhì)。而用虛擬組分來表征真實油品時，實際上等同于用圖中低密度的點去修正中心點，只是修正過程涉及了其區(qū)域中的全部點。

圖13 PC8中每個真實組分的局部密度與所在空間位置的可視化Fig.13 Display of local density and spatial location of each real component in PC8

2.6 精餾模擬驗證

下面使用ASPEN 進(jìn)一步驗證集總方法的性能。使用本文方法（DPC）和動態(tài)聚類法（DynC）分別得到的10個虛擬組分進(jìn)行常壓精餾模擬，并與真實組分的結(jié)果進(jìn)行對比。此外，還將動態(tài)聚類法得到的兩種不恰當(dāng)集總的結(jié)果也加入比較。進(jìn)料條件及塔參數(shù)設(shè)置如圖14所示，模擬結(jié)果的對比展示于圖15 中。顯然，使用虛擬組分替代真實油品進(jìn)行精餾模擬時，方法面臨了嚴(yán)峻挑戰(zhàn)：隨著D/F（塔頂采出量與進(jìn)料量之比）的增大，4 個集總結(jié)果的塔頂采出溫度均被高估。其中DynC（10 Comp）方法的塔頂溫度MAE 為5.03℃，塔釜溫度MAE 為4.34℃，相比之下DPC 方法取得了更可令人接受的結(jié)果，其塔頂溫度的MAE 為1.95℃，精度提高61.2%；塔釜溫度的MAE 為2.01℃，精度提高53.7%，這為下一步的“組分還原”奠定了相對可靠的基礎(chǔ)。

圖14 進(jìn)料條件及精餾塔參數(shù)設(shè)置Fig.14 Feeding conditions and distillation column parameter settings

圖15 不同集總結(jié)果下的精餾計算結(jié)果Fig.15 Calculation results of distillation of different lumped results

2.7 模擬計算效率

當(dāng)組分?jǐn)?shù)目增加時，精餾塔的待求方程組規(guī)模也隨之增大，故采用少數(shù)虛擬組分代替大量真實組分進(jìn)行模擬計算和調(diào)和設(shè)計時，在計算效率上無疑是有明顯進(jìn)步的。下面以ASPEN 為模擬計算工具，使用微型計算機(jī)對10 個虛擬組分和287 個真實組分的精餾計算時間(CPU 使用時間)進(jìn)行比較（針對第2.6節(jié)中的流程），計算機(jī)的硬件配置如下。

通過MATLAB 與ASPEN 的通信功能獲取到兩種情況的計算時間(含程序的打開/關(guān)閉時間，每次計算前均保證CPU 占用一致)如下：287 個真實組分為47.856 s，10 個虛擬組分為9.671 s。顯然兩者存在顯著差距。此外，若過程模擬涉及包含精餾塔的撕裂流收斂，或針對精餾塔進(jìn)行優(yōu)化，則計算效率的差距將更大。

3 結(jié) 論

將基于密度峰的聚類算法思想引入油品分子集總過程，實現(xiàn)了集總數(shù)目（中心組分）的自動偵測和脫硫汽油組成的可靠劃分，并以某裝置的脫硫汽油實例進(jìn)行驗證，得到以下結(jié)論。

（1）在油品分子集總過程中，存在一個兼顧效率和準(zhǔn)確性的集總數(shù)目，而本研究提出的方法借助決策圖即可簡便地確定這個集總數(shù)目和少數(shù)離群組分，這在一定程度上克服了估計集總數(shù)目時的主觀性和盲目性。此外，該方法可對油品性質(zhì)數(shù)據(jù)集中分布不規(guī)則的數(shù)據(jù)或“小簇”進(jìn)行有效挖掘，在一定程度上避免了傳統(tǒng)劃分方法造成的分配錯誤，從而使油品分子得到更精確的劃分。

（2）中心組分代替真實油品的能力是有限的；而虛擬組分可綜合其子組分（真實組分）的性質(zhì)，更好地反映真實情況。但中心組分作為真實組分，其熱力學(xué)計算相比于虛擬組分而言是嚴(yán)格的，關(guān)于中心組分描述真實油品的能力是否可以優(yōu)化，還有待進(jìn)一步討論。

（3）實例驗證結(jié)果表明，DPC 集總方法得到的虛擬組分可以滿足相平衡計算的要求，該方法對于將來構(gòu)建完整的石油組分分子水平劃分具有較大潛力。

本工作在集總性質(zhì)的選取方面，現(xiàn)階段使用了固定的熱力學(xué)性質(zhì)組合及性質(zhì)權(quán)重，并未考慮是否有更優(yōu)的選擇。考慮到不同煉油工藝過程可能對所需的性質(zhì)組合和權(quán)重值要求不同，故之后的工作重點在于，通過優(yōu)化性質(zhì)組合提高集總結(jié)果對油品的可描述性以及不同需求下的最優(yōu)權(quán)重設(shè)計。

符 號 說 明

c——虛擬組分中的子組分總數(shù)

dc——截斷距離

dij——數(shù)據(jù)點i與數(shù)據(jù)點j的相互距離

MAE——平均絕對誤差

MW——分子量

N——集總性質(zhì)個數(shù)

n——油品真實組分總數(shù)

Pb——正常沸點下的壓力，MPa

Pbr——正常沸點下對比飽和蒸氣壓

Pc——純組分臨界壓力，MPa

pik——數(shù)據(jù)點i的第k個維度（性質(zhì)）的值

R——氣體常數(shù)，8.314 J/(mol·K)

Sg——相對密度

Tb——沸點，K

Tbr——對比正常沸點

Tc——純組分臨界溫度，K

wk——維度（性質(zhì)）k的權(quán)重值

γi——ρi與δi的乘積

δi——第i個數(shù)據(jù)點的相對距離

ρi——第i個數(shù)據(jù)點的局部密度

ω——偏心因子

下角標(biāo)

i,j——組分的樣本序號

k——性質(zhì)的樣本序號