汽油精制過程中的辛烷值損失預(yù)測模型

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

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

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

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

1 變量降維

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

1.1 采用改進(jìn)的因子分析降維

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

因子分析的主要步驟為:

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

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

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

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

(1)

其中,xi為標(biāo)準(zhǔn)化后的原始變量,fi為因子變量,且k

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

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

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

1.1.1 提取公因子方差

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

表1 公因子方差

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

1.1.2 提取公因子

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

表2 改進(jìn)的因子分析解釋的總方差

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

圖1 改進(jìn)后的因子分析法碎石圖

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

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

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

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

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

1.2 采用隨機(jī)森林算法降維

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

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

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

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

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

4)重復(fù)2)、3),直到所有節(jié)點(diǎn)都訓(xùn)練過了或者被標(biāo)記為葉子節(jié)點(diǎn)。

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

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

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

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

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

1.2.2 變量選擇

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

圖3 輸入變量重要性測度指標(biāo)柱形圖

由柱狀圖可知,模型中較為重要的有13個(gè)變量(由于數(shù)據(jù)量較大,橫坐標(biāo)的變量名稱顯示不完整),這里給出通過自查得出的結(jié)果,這13個(gè)變量分別為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。

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

圖4 變量對預(yù)測錯(cuò)誤率的影響

從對變量對預(yù)測錯(cuò)誤率的影響的角度看,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個(gè)較為重要的變量和原料的辛烷值共14個(gè)變量作為模型的主要變量。將此次降維得到的結(jié)果命名為結(jié)果2。

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

2 建立預(yù)測模型

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

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

圖6 隨機(jī)森林基本思想

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

利用隨機(jī)森林的預(yù)測過程如下:

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

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

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

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

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

將辛烷值損失設(shè)為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編程預(yù)測樣本點(diǎn)301～325的RON損失值,然后對比樣本觀測值來驗(yàn)證預(yù)測結(jié)果的準(zhǔn)確性,結(jié)果如表5。

表5 隨機(jī)森林預(yù)測結(jié)果

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

圖7 預(yù)測值與真實(shí)值的對比

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

3 總 結(jié)

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

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