基于機器學習的汽油精制工藝操作方法優化

張棟, 林建新, 劉博, 林坤

(北京建筑大學北京市城市交通基礎設施建設工程技術研究中心， 北京 100044)

隨著車用汽油質量標準的不斷提高，迫切需要提高汽油清潔化工藝[1]。汽油清潔的重點是降低汽油中的硫和烯烴含量，同時盡可能保持其辛烷值(RON)含量[2-3]。因此，對汽油精制化的工藝操作方法進行優化具有重要意義。

汽油中的辛烷值是反映汽油燃燒性能的最重要指標[4]。秦強等[5]利用隨機森林算法進行汽油辛烷值的近紅外光譜建模，以此提高預測模型精度。但近紅外光譜法對儀器的要求較高，在企業中難以實現覆蓋。部分學者也從構建數據關聯模型出發進行辛烷值預測。Anderson等[6]將乙醇納入汽油辛烷值關聯模型之中，利用實驗擬合的方式預測辛烷值，但模型擴展性和機制解釋性較差。Westbrook等[7]采用正向預測法，在特定條件下對目標燃料依據其放熱規律確定的其中辛烷值，但存在計算成本高昂的缺點。馬寅杰等[8]在此基礎上利用辛烷值正向預測原理，采用多項式混沌展開(polynomial chaos expansion, PCE)來代替正向預測方案，研究初始溫度、壓力及汽油典型成分對辛烷值得影響，提高了預測精度。雖然當前對辛烷值的測量方法取得一定的成果，但傳統的數據關聯模型中具有變量相對較少、機理建模對原料的分析要求較高，對過程優化的響應不及時等缺點。隨著傳感器技術的使用，在脫硫過程中可以收集大量的設備運行數據[9-10]。以此為基礎，利用數據挖掘技術，選取適當的機器學習算法[11-12]，可從大量設備運營數據中提取影響辛烷值預測的主要特征，優化預測精度。李煒等[13]結合混合粒子群遺傳算法BP(back propagation)神經網絡，對辛烷值進行預測，結果表明預測性能更優。石翠翠等[14]通過偏最小二乘回歸(partial least square regression, PLS)和互信息(mutual information, MI)組合算法選取與辛烷值相關的主要變量，然后利用改進天牛須搜索算法(improved long-horned beetle antennae search algorithm, ILBAS)優化BP神經網絡模型，從而提高辛烷值預測準確度。雖然各學者利用機器學習的方法提高了辛烷值的預測精度，但汽油精制中硫含量對汽油品質也產生了很大影響。因此如何在精制過程中通過改變操作變量達到抑制硫含量的同時保持辛烷值損失較小的研究亟待深入。

在對各操作變量開展特征降維選取主要變量的基礎上，分別構建辛烷值損失預測模型與硫含量預測模型，以此建立以硫含量最低為約束，辛烷值損失最低為目標的工藝操作方法優化模型，并設計遺傳算法求解得到不同的優化值，揭示辛烷值與硫含量的分布規律，以達到最優工業操作條件。

1 數據采集與預處理

為開展汽油辛烷值損失模型構建研究工作、研究選取2017年4月—2019年9月和2019年10月—2020年5月兩個時間段，共約3年的某石化企業的催化裂化汽油精制脫硫裝置運行歷史數據。原始數據[15]包括325個數據樣本、7個原料性質、2個待生吸附劑性質、2個再生吸附劑性質，以及另外354個操作變量，包括氫油比、反應過濾器壓差、還原器壓力、還原器流化氫氣流量、反應器上部溫度等基礎操作變量，具體數量如表1所示。

原始數據在采集過程中由于設備裝置出現錯誤等原因導致出現部分變量的部分數據為空值、數據異常等問題，致使數據質量較低。因此，有必要進行數據清洗以消除低質量數據所帶來的誤差。

表1 數據預處理范圍

數據清理步驟如下：

(1)填充缺失數據。對實驗中樣本全部空值的樣本進行刪除，對部分數據為空值的點位采用前后數據的平均值進行替換，然后對數據異常值進行處理。

(2)噪聲處理。根據工藝要求與操作經驗，利用數據處理中常用四分位距(interquartile range, IQR)[16]檢測異常值。四分位距就是上四分位與下四分位之間的差值，通過四分位距的 1.5 倍距離為標準，規定：超過上四分位加上1.5倍IQR距離，或者低于下四分位距減去1.5 倍IQR距離的點為異常值。

(3)3σ準則去除異常值。3σ準則[17]：設對被測量變量進行等精度測量，得到x1,x2,…,xn，算出其算術平均值x及剩余誤差vi=xi-x(i=1,2,…,n)，并按貝塞爾公式算出標準誤差σ，若某個測量值xb的剩余誤差vb(1≤b≤n)，滿足|vb|=|xb-x|>3σ，則認為xb是含有粗大誤差值的壞值，應予剔除。貝塞爾公式為

(1)

2 辛烷值損失預測模型

2.1 模型主要變量篩選

由于樣本數據操作變量維度高，加之較多操作變量會使模型精度下降，變量信息產生重疊，增大問題復雜度。因此，去除冗余信息，把高位數據降維，能有效減少冗余特征導致的重疊信息。研究首先通過三種方法計算各個變量與產品辛烷值之間的相關性，然后通過三種方法求得的相關性結果，選出主要變量。

2.1.1 皮爾森(Pearson)相關性系數

Pearson相關系數[18]是用來反映兩個變量之間相似程度的統計量，可用于計算特征與類別間的相似度，判斷所提取到的特征和類別是正相關、負相關還是沒有相關程度。

Pearson相關系數定義為兩個變量之間協方差和兩者標準差乘積的比值，相關的計算表達式為

(2)

2.1.2 斯皮爾曼(Spearman)相關性系數

Pearson相關性系數不考慮數據取樣點之間的距離，會使較小范圍內存在聯系的兩個采樣數據被整體范圍內數據的非相關性所掩蓋，但斯皮爾曼相關性系數[19]可根據原始數據的排序位置進行求解，解除這一限制。

Spearman相關系數在對于變量之間的排序具有較高的要求，位置關系會影響相關性程度的判斷。計算公式為

(3)

其中，皮爾森、斯皮爾曼相關性系數表示獨立變量與相關變量之間相互聯系的程度與相關性方向。兩者系數的取值范圍為[-1,1]。當相關性系數的絕對值越接近1時，證明兩隨機變量之間的聯系越強，而相關性系數的絕對值越接近0時，證明兩隨機變量之間的聯系越弱；當相關性系數為正時，表示兩者之間正相關，即一方變量或者大小往積極方向變化時，另一方的數值也會往積極方向變化，相關性系數為負時則相反。

2.1.3 最大信息系數

最大信息系數(maximal information coefficient, MIC)是一類可用于分析評估變量間代數關系和依賴程度的相關性算法，具備良好的通用性質和穩定性[20]。該系數能夠分析組合變量數據的相關性和內在聯系，解釋隱藏的相關性信息。

最大信息系數將各個變量兩兩進行處理，將獲得的數據進行歸一化處理，即可得到各個變量之間的相關性程度。由于是計算兩變量的信息疊加程度，因此該系數的取值范圍為[0,1],系數越接近1，則說明兩者間聯系越強烈。

2.1.4 確定主要變量

通過皮爾森相關性系數、斯皮爾曼相關性系數與最大信息系數量化各個變量對產品辛烷值的相關性。通過python以及MATLAB數據庫中的minepy數據包以及上述相關的計算函數，對所有的變量兩兩進行了分析，研究各個變量和產品辛烷值的相關性，計算獲得所有變量的相關性系數，對其進行排序，挑選相關性較強的前30個變量(變量編號所表示的含義見文獻[15])進行對比，如圖1和圖2所示。

OL表示烯烴體積分數，SH表示飽和烴(烷烴+環烷烴)體積分數， SC表示硫含量，RON表示辛烷值圖1 辛烷值RON與變量間Pearson系數Fig.1 Octane number RON and Pearson coefficient

OL表示烯烴體積分數，SH表示飽和烴(烷烴+環烷烴)體積分數， SC表示硫含量，RON表示辛烷值圖2 辛烷值RON與變量間Spearman系數Fig.2 Octane number RON and Spearman coefficient

結果顯示這些變量在皮爾森相關性系數和斯皮爾曼相關性系數下的正負相關并未改變，但系數值稍有不同，證明分析數據具有較高的數據質量，可用于數據分析。對各個變量進一步分析發現存在21個變量在這兩個系數下都有較強的相關性，如表2所示。

考慮到斯皮爾曼相關系數僅能表示出單調函數的非線型關系，而變量與目標值之間可能存在更加復雜的非線性關系。因此計算該21個主要變量與目標值之間最大信息系數。由圖3所示，僅原材料中的辛烷值RON的量與目標值的最大信息系數較大，而且其他變量系數值較小且并沒有明顯不同，說明各個參數與目標值之間并未有明顯的周期性關系。綜上所述，所挑選的21個主要變量與產品辛烷值具有強相關性。

表2 主要變量及相關系數

由表2可知，在挑選的21個變量中，原料中辛烷值對于產品辛烷值影響程度最大，該結論與實際結果相符合。此外，圖4是21個主要變量各自之間的相關性關系圖，多數變量的相關性系數絕對值處于0.35～0.55，結果表明，除原料辛烷值與產品辛烷值的相關性較強外，其他主要變量對于產品辛烷值的影響均處于同一水平，按照等級程度分級，兩者處于中等程度相關的級別。

OL表示烯烴體積分數，SH表示飽和烴(烷烴+環烷烴)體積分數， SC表示硫含量，RON表示辛烷值圖4 主要變量之間相關性關系Fig.4 Correlation between main variables

2.2 辛烷值損失預測模型

根據2.1節中選擇的21個主要變量，通過對各類機器學習中的預測模型對比，最終選擇以XGBoost[21]為基礎構建回歸模型進行辛烷值損失預測。

XGBoost是在原有梯度算法的基礎上進行了改進。模型為樹集成模型，通過將各個決策樹的決策結果進行總和，作為最終的預測值，即

(4)

梯度提升決策樹(gradient boosting decision tree, GBDT)通過上一次訓練得到的預測值殘差作為下一次訓練的標準值，能夠有效逼近樣本真實值，其主要原因在于GBDT使用樣本平方差損失作為損失函數，其基本表達式為

(5)

而XGBoost則在GBDT基礎上將損失函數從泰勒展開的第一階設定為第二階，使用前兩階作為改進的殘差。

此外，為了限制模型的復雜程度，使決策樹模型的各個決策樹處于“弱勢”狀態，保證最終結果不受異常訓練集的影響，加入正則化項，XGBoost使用葉子節點的個數作為替代指標來減低模型的復雜程度，均衡各個決策樹預測結果的權重。正則化項公式為

(6)

式中：λ、γ表示人工設置參數；T為葉子節點數；w為決策樹所有葉子節點值形成的向量。

因此，其目標函數由梯度提升算法損失和正則化項構成，公式為

(7)

2.3 交叉驗證

針對樣本數據不足、各個變量之間差異較大以及維度高，研究選取交叉驗證[22]的方法進行分類訓練，如圖5所示?；静襟E如下。

(1)將數據樣本均分為10等份。

(2)選取樣本1作為測試集數據，樣本2～10作為訓練集數據進行第一次回歸預測。

(3)選取樣本2作為測試機數據，樣本3～10作為訓練集數據進行第二次回歸預測。

(4)依此類推，通過不同的測試集數據與訓練集數據，獲得各次訓練的擬合優度，對其求平均值，即當作該操作變量的回歸預測的結果。

(5) 使用交叉驗證方法對所有操作變量進行回歸預測，進行分析。

圖5 交叉驗證示意圖Fig.5 Schematic diagram of cross-validation

3 工藝操作方法優化

為提高汽油精制化工藝，在產品硫含量不大于5 μg/g的前提下，研究利用建立的損失預測模型，對獲得辛烷值損失降幅大于30%的樣本的主要變量進行優化。

首先，建立硫含量預測模型：研究假設與辛烷值損失量相關性較強的操作變量對硫含量亦有較強相關性，同時在以硫含量是否超過5 μg/g為分界線建立硫含量預測模型，同時驗證預測模型有效性。其次，在保證硫含量滿足要求的情況下，以2.1節中的主要變量作為決策變量，利用遺傳算法-聚類遞歸的方法對辛烷值損失降幅進行優化，并提取辛烷值損失大于30%的樣本下的最優操作條件。具體流程如圖6所示。

圖6 優化流程圖Fig.6 Optimization flow chart

3.1 硫含量預測模型

根據高維度樣本量較少且選取變量相關性較強的特點，研究選取嶺回歸[23]作為硫含量預測模型。嶺回歸是一種專用于共線性數據分析的有偏估計回歸方法，其使一種改良的最小二乘估計法，通過放棄最小二乘法的無偏性，來解決高維度樣本量的“病態”數據擬合問題。

嶺回歸是在線性回歸的基礎上加入了正則項，增加了模型的泛化能力。正則項一般采用一、二范數，并解決線性回歸總不可逆的情況。表達式為

(8)

式(8)中：λ為嶺系數；θ為回歸系數；hθ(xi)為預測值。

3.2 遺傳算法尋優

為通過預測模型獲得使辛烷值損失量最小的操作變量條件，選取遺傳算法作為優化算法對實驗中操作變量值進行求解。將每一組操作變量的組合稱作1條染色體，每個操作變量視為染色體的基因，在原料性質、待生吸附劑和再生吸附劑的性質數據保持不變的條件下，以辛烷損失最小作為目標函數。遺傳算法具體操作變量尋優過程如下所示：

(1)生成初代染色體：在354個操作變量信息[15]中各個操作變量的取值范圍之內，按照隨機均勻分為的原則隨機生成100個操作變量，即生成100個基因；并且相互結合形成100組操作變量組合，即100條染色體。

(2)將操作變量組合分別輸入辛烷值回歸預測模型以及硫含量預測模型，分別得出辛烷值預測損失量以及硫含量是否超過5 μg/g。

(3)選取實驗結果：計算辛烷值預測損失量，并觀察其變化是否平穩，如果連續2代辛烷值損失量之間預測差值穩定在0.01之內，則停止迭代，否則繼續進行下一步。

(4)基因突變與基因重組：按照辛烷值損失預測值由小到大排序，選擇辛烷值預測損失量的前20條染色體進行突變與重組。下一代染色體中突變染色體占40條，重組染色體占30條。突變規則為按照各個操作變量的Δ值(各主要操作變量每次允許調整幅度)幅度變化；重組規則為參考上一代中的前20條染色體基因庫，隨機選取一個值作為本代基因。

(5)轉向(2)重新開始實驗。

3.3 聚類遞歸

由于操作變量在優化的過程中逐漸聚類，不同操作變量聚類分布特性不同。遺傳算法結果顯示操作變量取值存在多組最優結果，如果取使辛烷值損失最小的操作變量取值組合，看似結果最優，但操作變量間關系復雜，若某一操作變量稍微改變就引起辛烷值損失較大，說明模型魯棒性不強。如果直接取均值，從各個操作變量的尋優結果分布特性來看，無法有效提升模型魯棒性，缺乏科學依據。為防止這一情況發生，研究提出聚類遞歸[24]開展操作變量取值，增強模型可移植性。

由于K-means算法需自行選取初始聚類中心且易陷入局部最優，因此研究選取DBSCAN算法通過密度可達性的性質將樣本數據集合分為K個簇標記數組以及數組外的噪聲點集合，公式為

(9)

對操作變量的取值集中程度進行量化評價，研究選取“分散度作為量化評價指標以計算操作變量的數據分布集中程度。對于單水平集中的操作變量具有較強的集中性。對于多水平集中的變量應先聚類，再通過量化評價指標計算其集中程度，集中性越強的操作變量可以優先取值。對于多水平集中的操作變量數據內部分類而言，可以首先選取占比較大的子類作為取值集合，計算各個子類的極值，以數據各子類的占比作為權重系數進行加和計算來表征數據整體的集中程度，稱之為分散度。具體計算過程如下：

(1)對各個操作變量以操作變量信息中操作范圍為準對結果數據進行歸一化處理，并進行聚類分析(DBSCAN)；

(2)計算每類數據的極值(δ)與占比(ρ)；

(3)按照公式計算各參數標定結果的分散程度，即

(10)

分散度小說明數據集中程度更強，按照分散度由小到大順序依次確定操作變量取值。對分散度最小的操作變量進行DBSCAN聚類，保留占比最大的子類，計算子類均值作為該操作變量的最終取值，依次確定各個操作變量。

4 優化方法應用分析

在構建辛烷值損失預測模型中，實驗將294個樣本數據作為訓練集、30個樣本數據作為測試集。通過網格搜索法，以模型預測測試集的均方誤差最小作為尋優目標，搜索參數最優組合。依托Python平臺，調取XGBoost庫，以2.1節選取的21個主要變量的樣本數據作為輸入，辛烷值損失值作為輸出，建立XGBoost回歸預測模型。通過Sklearn機器學習庫中Model_Selection模塊，使用GridSearchCV進行網格搜索。預測結果如表3所示。

表3 辛烷值模型預測結果

通過交叉驗證驗證模型有效性，模型最終交叉驗證平均準確率96.54%，R2擬合優度平均值為0.784。

為了得到硫含量的預測值，將樣本數據集分為30%測試數據、70%訓練數據，調用Python中Sklearn庫函數Ridge，對硫含量進行回歸分析。預測結果如圖7所示，預測值的平均絕對誤差為0.15，均方差為0.39，表明預測效果良好，可作為硫含量的預測模型。

隨著遺傳算法的不斷迭代，辛烷值損失值越來越趨于平穩。325個樣本中平均經過12.68次迭代，平均1 268次優化實驗，如圖8所示。

通過3.3節中聚類遞歸計算得到各個變量的分散度，計算結果如表4所示。以2號樣本為例，操作變量No.7(D104溫度)在迭代過程中逐漸集中，由圖9(a)可知，變量值更傾向于單一子類(黃色)，集中性較強。變量取值分布更接近于取值范圍(100～150 ℃)的上限，樣本中取值為125 ℃，最終優化結果142 ℃。反觀變量No.48(3#催化汽油進裝置流量)見圖9(b)，取值范圍為0～90 t/h，向兩個子類集中，選取占比較大的子類并取均值作為該參數最終取值，可知應為綠色子類均值，最終樣本中優化結果61 t/h。而圖9(c)、圖9(d)兩個操作變量，取值分布向多個子類集中，同樣選取占比最大子類取均值作為變量最終取值。

圖7 嶺回歸預測結果Fig.7 Ridge regression prediction results

圖8 遺傳算法迭代圖Fig.8 Iterative graph of genetic algorithm

表4 參數標定結果的分散度計算值

圖9 部分操作變量聚類遞歸分布Fig.9 Clustering recursive distribution of some operating variables

根據第3節工藝操作要求，依次對辛烷值損失量優化降幅超過30%的308組樣本進行操作變量的聚類遞歸，最終優化方案平均降幅達50.26%。這充分說明了本文提出優化流程的可靠性。

以樣本1為例，優化前后的各個參數取值對比如圖10所示?？梢钥闯觯饕獌灮瘲l件為降低No.186(原料汽油硫含量)、No.26(輕烴出裝置流量)、No.48(3#催化汽油進裝置流量)。具體操作條件如表5所示。

在實際工業流程中，工業裝置為了平穩生產，優化后的主要操作變量往往只能逐步調整到位。研究以133號樣本為例，在原料性質、待生吸附劑和再生吸附劑的性質數據保持不變的情況下，調整主要變量，設每次允許調整的幅度值為Δ。為最快達到最優，各個變量取各自范圍內最大限度逐步逼近最優參數值，當到達最優值時，該變量取值不應再變化。通過數據分析計算，No.26(輕烴出裝置流量)到達最優值需要移動步點最多，為53步，No.55(原料換熱器管程總管進口溫度)需要移動步點最少，為2步；各主要變量尋優過程信息如表6所示。

圖10 樣本1優化前后對比圖Fig.10 Comparison diagram of sample 1 before and after optimization

將表6中各主要變量的數據值變化輸入已構建的辛烷值損失回歸預測模型與硫含量回歸預測模型來預測每一步辛烷值和硫含量的變化軌跡并生成可視化圖形，如圖11所示。

從圖11中可以看出，各個變量在逐步調節過程中，硫含量起伏變化相對較大，在第24步時雖然辛烷值損失不是最少，但硫含量已經降到最低，企業可根據需求選取適當的主要變量取值。

表6 變量優化過程信息表Table 6 Variable optimization process information table

圖11 辛烷值與硫含量變化軌跡Fig.11 Change track of octane number and sulfur content

5 結論

(1) 研究分別基于XGboost、嶺回歸構建辛烷值損失預測模型、硫含量預測模型，證明了預測模型的有效性與可行性。

(2) 形成采集數據預處理-樣本特征降維-預測模型構建-交叉驗證-遺傳算法求解-主要變量聚類回歸的工藝操作變量取值優化流程。

(3) 研究將不同操作變量進行DBSCAN聚類遞歸，計算多水平集下各操作變量的分散度，選取較大分散度的子類作為操作條件，大幅提高了各操作條件的取值準確性，使優化模型更加精準。

(4) 構建工藝操作取值可視化展示更加直觀反映了各主要變量的變化范圍以及對應變量下硫含量與辛烷值含量的變化趨勢。這為企業在實際生產中提供了更優的操作方案。