杜明洋,張?zhí)鹛?薄其高,許文文

齊魯工業(yè)大學(山東省科學院) 數(shù)學與統(tǒng)計學院,濟南 250353

目前,小型車輛的主要燃料是汽油,汽油燃燒排放的產(chǎn)物主要有一氧化碳、氮氧化合物、碳氫化合物、微粒等,對大氣環(huán)境造成了重大影響。為此,世界各國都制定了日益嚴格的汽油質(zhì)量標準。汽油清潔化重點是降低汽油中的硫、烯烴含量,同時盡量保持其辛烷值。為滿足日益嚴格的清潔汽油標準不斷降低硫和烯烴含量的需求,國內(nèi)外在汽油清潔化領域開展了大量的研究工作。

辛烷值(RON)是反映汽油燃燒性能的最重要指標,目前除考慮改進發(fā)動機的構(gòu)造和操作條件外,主要依靠提高汽油辛烷值,借以提高發(fā)動機壓縮比的辦法來增加功率和節(jié)約汽油。但是現(xiàn)有技術在對催化裂化汽油進行脫硫和降烯烴過程中,普遍降低了汽油辛烷值。辛烷值每降低1個單位,相當于損失約150元/噸。以一個100萬噸/年催化裂化汽油精制裝置為例,若能降低0.3個單位RON損失,其經(jīng)濟效益將達到4 500萬元。因此降低汽油辛烷值損失具有重要意義。本文采用改進的因子分析法和隨機森林算法建立汽油精制過程中的辛烷值損失預測模型。具體數(shù)據(jù)來源于某石化企業(yè)的催化裂化汽油精制脫硫裝置運行4年所積累的歷史數(shù)據(jù)[1-3]。

1 變量降維

在解決實際問題時,自變量過多對于模型的構(gòu)建是弊大于利的,因此,一般采用先降維再建立模型的方法來簡化模型。本文采用改進的因子分析法、隨機森林兩種方法分別對操作變量進行降維,最后比較兩種方法的降維效果。

1.1 采用改進的因子分析降維

因子分析是將多個實測變量轉(zhuǎn)換為少數(shù)幾個綜合指標(或稱因子變量),它反映一種降維的思想。通過降維將相關性高的變量聚在一起,從而減少需要分析的變量的數(shù)量,降低問題分析的復雜性。因子變量滿足以下特點:個數(shù)遠遠少于原始變量,但能反映原始變量的絕大部分信息。

因子分析的主要步驟為:

1)判斷待分析因子變量是否適合進行因子分析。

2)構(gòu)造因子變量并計算公因子方差。

3)利用旋轉(zhuǎn)法使因子變量具有可解釋性。

因子分析的數(shù)學模型為:

(1)

其中,xi為標準化后的原始變量,fi為因子變量,且k

上述模型也可用矩陣表示,即A=AF+ε,其中F為因子變量,A為載荷矩陣,aij為因子載荷,ε為特殊因子。所謂因子載荷,是指在因子變量不相關的情況下,第i個原始變量與第j個原始變量的相關系數(shù),aij的絕對值越大,則xi與fi的相關性越強[4-5]。

將包括7個原料性質(zhì)、2個待生吸附劑性質(zhì)、2個再生吸附劑性質(zhì)、2個產(chǎn)品性質(zhì)以及另外354個操作變量(共367個變量)作為變量,RON損失作為因變量進行因子分析。將因變量和操作變量進行編號,令辛烷值損失為Y,367個變量編號依次為x1,x2,…,x367。

由于工業(yè)數(shù)據(jù)相對比較復雜,本文在傳統(tǒng)因子分析的基礎上,對原始操作變量與因變量進行相關性分析,分別計算原始操作變量與因變量之間的相關性,提取出與因變量相關性較強的變量,然后進行因子分析。此方法可以減少因子分析的工作量,同時提高降維效果,利用SPSS對原始數(shù)據(jù)進行分析處理。

1.1.1 提取公因子方差

每個操作變量均能通過公因子表達,但每個公因子可以表達原始數(shù)據(jù)信息的多少是通過公因子方差表的“提取”來表示的。利用SPSS得到公因子方差表,如表1。

表1 公因子方差

由上表可知,所有變量的提取值均大于0.7,說明這些變量都可以很好的被公因子表達,即分析結(jié)果比較理想,原始變量適合進行因子分析。

1.1.2 提取公因子

對預處理后的數(shù)據(jù)中原始操作變量與因變量進行相關性分析,得到相關性矩陣,將相關系數(shù)0.1作為標準,剔除相關性低于0.1的變量,共剔除133個變量。對保留的234個變量的數(shù)據(jù)進行因子分析,利用SPSS輸出得到解釋的總方差表如表2。

表2 改進的因子分析解釋的總方差

將特征值大于1因子作為公因子,由表2可知,可提取14個公因子,其累計方差貢獻率為84.779%,即這14個公因子可以解釋84.779%的變量信息,原有變量的大多數(shù)信息能很好的被這14個公因子涵蓋。輸出的碎石圖如圖1。

圖1 改進后的因子分析法碎石圖

左側(cè)斜率較大的部分代表特征值大,即成份1、成份2、…、成份14這14個成份包含大部分變量信息,提取這14個成份為公因子是可行的。

通過方差最大法獲得正交旋轉(zhuǎn)因子載荷矩陣,如表3。

表3 改進的因子分析旋轉(zhuǎn)成分載荷矩陣

根據(jù)旋轉(zhuǎn)成分矩陣可知,第1個公因子上具有較大載荷值的是原料性質(zhì)中的辛烷值、S-ZORB.FT_9301.TOTAL、S-ZORB.FT_9201.TOTAL、S-ZORB.FT_9403.TOTAL和S-ZORB.FT_3301.TOTAL；第2個公因子上具有較大載荷值的是S-ZORB.FT_1503.TOTALIZERA.PV和S-ZORB.FT_1504.TOTALIZERA.PV；第3個公因子上具有較大載荷值的是S-ZORB.PT_5201.DACA；第4個公因子上具有較大載荷值的是S-ZORB.PDT_2703B.DACA和S-ZORB.PDI_2801.DACA；第5個公因子上具有較大載荷值的是S-ZORB.FT_1001.PV；第6個公因子上具有較大載荷值的是S-ZORB.FT_9102.PV；第7個公因子上具有較大載荷值的是S-ZORB.PT_2901.DACA；第8個公因子上具有較大載荷值的是S-ZORB.TE_7106B.DACA；第9個公因子上具有較大載荷值的是S-ZORB.TE_1102.DACA.PV；第10個公因子上具有較大載荷值的是S-ZORB.PDC_2702.DACA；第11個公因子上具有較大載荷值的是S-ZORB.PT_6005.DACA；第12個公因子上具有較大載荷值的是S-ZORB.FC_5103.DACA；第13個公因子上具有較大載荷值的是S-ZORB.FT_9001.PV；第14個公因子上具有較大載荷值的是S-ZORB.TC_3102.DACA。

因此,根據(jù)改進的因子分析法進行降維提取出的公因子的信息大致可以用這20個變量包含,將此降維結(jié)果命名為結(jié)果1。

1.2 采用隨機森林算法降維

隨機森林是隨機的建立一個森林,森林由很多的決策樹組成,隨機森林的每一棵決策樹之間沒有關聯(lián)。得到森林后,當有一個新的輸入樣本進入時,讓森林中的每一棵決策樹分別判斷這個樣本應該屬于哪一類(對于分類算法),然后觀察哪一類被選擇最多,就預測這個樣本為該類。隨機森林的訓練過程總結(jié)如下:

1)給定訓練集S,測試集T,特征維數(shù)F。確定參數(shù):決策樹的數(shù)量t,每棵樹的深度d,每個節(jié)點使用到的特征數(shù)量f。終止條件:節(jié)點上最少樣本數(shù)s,節(jié)點上最少的信息增益m。

對于第1-t棵樹,i=1-t:

2)從S中有放回的抽取大小和測試集一樣的訓練集S(i),作為根節(jié)點的樣本,從根節(jié)點開始訓練。

3)如果當前節(jié)點上達到終止條件,則設置當前節(jié)點為葉子節(jié)點；如果是分類問題,該葉子節(jié)點的預測輸出為當前節(jié)點樣本集合中數(shù)量最多的那一類c(j),概率m為c(j)占當前樣本集的比例；如果是回歸問題,預測輸出為當前節(jié)點樣本集各個樣本值的平均值。然后繼續(xù)訓練其他節(jié)點。如果當前節(jié)點沒有達到終止條件,則從F維特征中無放回的隨機選取f維特征。利用這f維特征,尋找分類效果最好的一維特征k及其閾值th,當前節(jié)點上樣本第k維特征小于th的樣本被劃分到左節(jié)點,其余的被劃分到右節(jié)點。繼續(xù)訓練其他節(jié)點。

4)重復2)、3),直到所有節(jié)點都訓練過了或者被標記為葉子節(jié)點。

5)重復2)、3)、4),直到所有CART都被訓練過[6-8]。

此處運用R軟件對數(shù)據(jù)通過隨機森林的方法進行處理。

1.2.1 確定決策樹數(shù)量

首先尋找最佳參數(shù)ntree,即指定隨機森林所包含的最佳決策樹數(shù)目。通過R軟件編程輸出模型誤差與決策樹數(shù)量關系圖如圖2,易見當決策樹數(shù)量大于等于300時,模型誤差趨于穩(wěn)定。因此選擇300作為最佳決策樹數(shù)量[9]。

圖2 模型誤差與決策樹數(shù)量關系圖

1.2.2 變量選擇

為評價輸入變量的重要性,可用importance成分繪制對各輸入變量添加辛烷值后對整體預測精度平均影響的柱狀圖,如圖3。

圖3 輸入變量重要性測度指標柱形圖

由柱狀圖可知,模型中較為重要的有13個變量(由于數(shù)據(jù)量較大,橫坐標的變量名稱顯示不完整),這里給出通過自查得出的結(jié)果,這13個變量分別為S-ZORB.LC_5001.PV、S-ZORB.TE_5202.PV、S-ZORB.FC_5202.PV、S-ZORB.FT_9202.PV、S-ZORB.TE_1105.PV、S-ZORB.PT_6002.PV、S-ZORB.FT_9101.TOTAL、S-ZORB.FT_3201.DACA、S-ZORB.DT_2001.DACA、S-ZORB.PDT_2001.DACA、S-ZORB.TE_7506.DACA、S-ZORB.TE_2104.DACA.PV、S-ZORB.PC_1001A.PV。

為更全面的評價各輸入變量的重要性,可采用varImpPlot函數(shù)得到變量對預測錯誤率的影響排名及變量對GINI系數(shù)下降率的影響排名。變量對預測錯誤率的影響表示隨機去掉一個特征,模型準確度下降的百分比,越高表明特征越重要；變量對GINI系數(shù)下降率的影響表示使用某一個特征進行分裂時,GINI系數(shù)下降的平均幅度,越高表明該特征的分裂質(zhì)量越好,即該特征越重要。具體排名如圖4、圖5[10]。

圖4 變量對預測錯誤率的影響

從對變量對預測錯誤率的影響的角度看,S-ZORB.PC_1001A.PV、S.ZOBB.PDT_3502.DACA、S-ZORB.DT_2001.DACA、S-ZORB.TE_1107.DACA、S-ZORB.TE_7506.DACA、S-ZORB.TE_5202.PV、 S-ZORB.TE_1001A.PV、 S-ZORB.FT_9302.PV、S-ZORB.AT.0009.DACA.PV、 S-ZORB.TE_6008.DACA等較為重要。

圖5 變量對GINI系數(shù)下降率的影響

從變量對GINI系數(shù)下降率影響的角度看,辛烷值RON、 S-ZORB.TE_1001A.PV、S-ZORB.DT_2001.DACA、S-ZORB.TE_7506.DACA、S-ZORB.TE_5002.DACA、S-ZORB.TE_1201.PV、S-ZORB.LC_5101.PV、S.ZOBB.PDT_3502.DACA、性質(zhì)辛烷值RON、S.ZOBB.PDP_35021003.DACA等變量對辛烷值損失的影響較大。

綜合考慮三種角度提取的主要變量,選擇S-ZORB.LC_5001.PV、S-ZORB.TE_5202.PV、S-ZORB.FC_5202.PV、S-ZORB.FT_9202.PV、S-ZORB.TE_1105.PV、S-ZORB.PT_6002.PV、S-ZORB.FT_9101.TOTAL、S-ZORB.FT_3201.DACA、S-ZORB.DT_2001.DACA、S-ZORB.PDT_2001.DACA、S-ZORB.TE_7506.DACA、S-ZORB.TE_2104.DACA.PV、S-ZORB.PC_1001A.PV將這13個較為重要的變量和原料的辛烷值共14個變量作為模型的主要變量。將此次降維得到的結(jié)果命名為結(jié)果2。

對比兩種降維方法所得結(jié)果的擬合優(yōu)度R2,結(jié)果2的R2較大,即隨機森林的降維效果較好。

2 建立預測模型

沿用隨機森林篩選出的主要建模變量,利用R語言編程分析降維后的變量與RON損失的線性關系,建立RON損失預測模型,并預測最后25個時間點的RON損失值。最后,通過預測結(jié)果與實際值的對比驗證RON損失預測模型的可靠性[11]。

隨機森林回歸模型的確立主要分為兩個步驟,首先是決策樹的構(gòu)建,其次是森林的形成。隨機森林回歸是由很多回歸決策樹模型{h(x,θk),k=1}組成的組合分類模型,且參數(shù)θk是獨立分布的隨機向量,在給定自變量X下,每個決策樹回歸都會有一個預測結(jié)果。它的基本思想與流程如圖6。

圖6 隨機森林基本思想

其中，X:初始樣本,Xi*:重抽樣樣本。

利用隨機森林的預測過程如下:

對于第S棵樹,i=1-t:

1)從當前樹的根節(jié)點開始,根據(jù)當前節(jié)點的閾值th,判斷是進入左節(jié)點

2)重復執(zhí)行1)直到所有樹都輸出了預測值。如果是分類問題,則輸出為所有樹中預測概率總和最大的那一個類,即對每個c(j)的p進行累計；如果是回歸問題,則輸出為所有樹的輸出的平均值。

針對本文的研究對象,將前300個時間位點的數(shù)據(jù)作為訓練數(shù)據(jù),第301～325時間位點的數(shù)據(jù)作為檢測數(shù)據(jù)即進行驗證,因此先將降維后的主要建模變量利用R軟件進行線性回歸,輸出結(jié)果如表4。

表4 隨機森林回歸結(jié)果

將辛烷值損失設為y,根據(jù)表中數(shù)據(jù)可得回歸方程如下:

y=-8.928 0+0.006x2+0.008 2x30+

0.021 5x33+0.001 5x34+0.000 2x49-

0.000 3x68-0.425 8x78-0.000 3x191+0.001 9x246-0.006x249+0.019 1x338-

0.004 6x367,

(2)

根據(jù)上述回歸方程通過R編程預測樣本點301～325的RON損失值,然后對比樣本觀測值來驗證預測結(jié)果的準確性,結(jié)果如表5。

表5 隨機森林預測結(jié)果

再根據(jù)表中數(shù)據(jù)通過Excel進行作圖,如圖7。

圖7 預測值與真實值的對比

從表5和圖7可以看出,預測值與真實值的差別不大,因此可以認定本文所得RON損失預測模型是可靠的。

3 總 結(jié)

本文研究了一般石化企業(yè)的催化裂化汽油精制脫硫裝置的模型構(gòu)建。首先基于處理后的數(shù)據(jù)綜合運用改進的因子分析法和隨機森林法兩種方法通過SPSS、R軟件分別對變量進行降維處理并對比降維效果。其次,通過R語言編程分析降維后的變量與RON損失的線性關系,建立RON損失預測模型,并預測最后25個樣本點的RON損失值。最后,通過預測結(jié)果與實際值的對比驗證了RON損失預測模型的可靠性。

該問題的研究具有重要的理論意義和實際應用價值,后續(xù)還可進一步尋求其他模型檢驗以期達到更好的效果。