陸天浩，李玲玲，陳宇峰，湯 瓊

0 引言

汽油的辛烷值作為汽油的商品牌號（例如89#、92#、95#），是衡量汽油發(fā)動機(jī)燃料抗爆性能優(yōu)劣的重要指標(biāo)。辛烷值越高，表示其抗爆性能越好，發(fā)動機(jī)壓縮比越高[1]。辛烷值的提高可以提升發(fā)動機(jī)功率，增加車輛行程；同時還能起到節(jié)約燃料，減輕使用者經(jīng)濟(jì)負(fù)擔(dān)的作用。對于小型車輛而言，汽油是其主要的燃料，汽油燃燒產(chǎn)生的尾氣是影響大氣環(huán)境的重要因素。辛烷值反映了汽油的燃燒性能，然而當(dāng)前的脫硫催化裂化汽油技術(shù)，使汽油的辛烷值大大降低。因此提高辛烷值是提高經(jīng)濟(jì)效益和改善環(huán)境的一個重要手段[2]。

為了解決上述問題，針對“華為杯”第十七屆中國研究生數(shù)學(xué)建模競賽①本文為“華為杯”第十七屆中國研究生數(shù)學(xué)建模競賽二等獎?wù)撐?。B 題：汽油辛烷值建模，將BP（back propagation）神經(jīng)網(wǎng)絡(luò)模型與遺傳算法[3-8]引入預(yù)測與優(yōu)化任務(wù)中，并通過預(yù)測優(yōu)化找到符合條件的優(yōu)化方案。

1 數(shù)據(jù)處理

1.1 數(shù)據(jù)預(yù)處理

用數(shù)學(xué)建模競賽B 題所提供的數(shù)據(jù)集作為原始數(shù)據(jù)集。在數(shù)據(jù)集中，共有325 個樣本，根據(jù)樣本采集時間排序，采集時間為2017-04—2020-05。每個樣本包含由原料性質(zhì)、產(chǎn)品性質(zhì)、待生吸附劑性質(zhì) 、再生吸附劑性質(zhì)等不同操作變量構(gòu)成的354 個采集位節(jié)點。由于是原始數(shù)據(jù)，不同樣本中的不同位節(jié)點存在未采集數(shù)據(jù)或數(shù)據(jù)缺失等情況。如果不同樣本的同一位節(jié)點缺失數(shù)據(jù)較多，則該位節(jié)點不能準(zhǔn)確反映變化情況，需要刪除該位節(jié)點。如果不同位節(jié)點的同一樣本數(shù)據(jù)較多，同樣需要刪除該樣本。對于缺失值個數(shù)較少的位節(jié)點，可根據(jù)臨近時間的樣本對應(yīng)位節(jié)點值進(jìn)行填充。再對篩選后的數(shù)據(jù)采用3σ準(zhǔn)則去除異常值；根據(jù)實際應(yīng)用背景確定取值范圍，將超過范圍的位節(jié)點或?qū)?yīng)樣本數(shù)據(jù)剔除。具體數(shù)據(jù)處理如下。

1）遍歷整個數(shù)據(jù)表，將采集值缺失過多的樣本及位節(jié)點刪除。經(jīng)過篩選，共刪除12 個位節(jié)點，具體名稱如表1 所示。這類位節(jié)點由于信息缺失過多，因此不能清晰體現(xiàn)該位節(jié)點對辛烷值損失的影響。

表1 第一步被刪除的位節(jié)點Table 1 Removed sites in the first step

2）對位節(jié)點進(jìn)行過濾操作，刪除數(shù)值全部為空值或無法通過相鄰樣本進(jìn)行填充的位節(jié)點。前者沒有數(shù)值變化，不能體現(xiàn)該位節(jié)點對辛烷值損失的影響。后者只有部分時間段樣本，不能完整地體現(xiàn)該位節(jié)點對辛烷值損失的作用。過濾后，刪除19 個位節(jié)點，具體名稱如表2 所示。

表2 第二步被刪除的位節(jié)點Table 2 Removed sites in the second step

3）對于存在少量空缺值的位節(jié)點，用相鄰樣本的均值進(jìn)行填補(bǔ)，需要補(bǔ)充的位節(jié)點如表3 所示。

表3 缺省值填補(bǔ)位節(jié)點Table 3 Nodes filled by the default value

4）根據(jù)不同操作變量實際取值范圍，剔除位節(jié)點采集值異常的32, 29, 27 號樣本。

5）使用3σ準(zhǔn)則去除異常值及其所帶來的影響。

完成上述預(yù)處理后，將285, 313 號樣本結(jié)果加入到對應(yīng)位置。

1.2 數(shù)據(jù)降維

模型1 主成分分析

主成分分析法是使用線性降維的方法來尋找模型的主要操作變量，以達(dá)到精簡影響因素數(shù)量的目的。本文使用SPSS 進(jìn)行主成分分析，使多維的數(shù)據(jù)降維得到少數(shù)幾個綜合指標(biāo)，具體操作步驟如下：

1）數(shù)據(jù)標(biāo)準(zhǔn)化

設(shè)x1,x2,…,xm為m個主成分分析變量，aij為m個評價對象中第i個評價對象對應(yīng)于第j個分析變量的取值。將每個aij值標(biāo)準(zhǔn)化，

標(biāo)準(zhǔn)化指標(biāo)變量，即

2）計算數(shù)據(jù)的相關(guān)系數(shù)矩陣

相關(guān)系數(shù)矩陣R的具體表達(dá)式為

式中，rpq為第p個指標(biāo)與第q個指標(biāo)之間的相關(guān)系數(shù)，且

使用主成分分析法對數(shù)據(jù)集中的320 個位節(jié)點進(jìn)行分析。第一次實驗中，將主成分個數(shù)設(shè)置為25，共篩選出200 多個變量。但由于決定主成分因子的操作變量種類各不相同，在工藝和設(shè)備操作上沒有很強(qiáng)的關(guān)聯(lián)性。因此將主成分個數(shù)設(shè)為5，然后重新實驗，結(jié)果如下：反應(yīng)過濾器壓差，反應(yīng)器上部溫度，SZORB.FT_9001.TOTAL，SZORB.FT_5201.TOTAL，SZORB.FT_9102.TOTAL。第二次實驗結(jié)果優(yōu)于第一次，但由于傳統(tǒng)方法的關(guān)聯(lián)模型變量少、要求高，因此響應(yīng)慢、效果差，且5 個主成分因子只能解釋59.374%的結(jié)果。

總的來說，通過對位節(jié)點進(jìn)行降維，可以得到一部分主要操作變量。但是由于操作變量之間具有高度非線性和相互強(qiáng)耦合的關(guān)系，使用主成分分析法篩選出的主要操作變量影響因子較低，不具有很強(qiáng)的解釋性，因此本文采用非線性的降維方法。

模型2 層次聚類法

由于數(shù)據(jù)集中的320 個有效位節(jié)點之間具有高度非線性和相互強(qiáng)耦合關(guān)系，為了尋找具有代表性和獨立性的主要操作變量，再采用非線性的層次聚類法進(jìn)行篩選。最后將主成分分析與層次聚類法篩選出的結(jié)果進(jìn)行綜合分析，得出最終建模需要的主要操作變量。構(gòu)建過程如下：

參考?xì)W幾里德距離評價，在數(shù)據(jù)不規(guī)范時，皮爾遜相關(guān)系數(shù)會給出更好的結(jié)果，故本文使用皮爾遜相關(guān)系數(shù)來尋找模型的主要操作變量。具體算法見式（5）（6）。

式（5）（6）中：X、Y為兩個不同的變量；

μX、μY和D(X)、D(Y)分別為變量X、Y的均值和方差；

ρXY為兩變量X、Y間的皮爾遜相關(guān)系數(shù)，當(dāng)2個變量完全匹配時，ρXY=1，當(dāng)2 個變量毫無關(guān)系時，ρXY=0。

為了使相似度越大的兩個變量之間距離越小，采用1 與皮爾遜相關(guān)系數(shù)的差值來衡量[9]，見式（6）。

使用SPSS 進(jìn)行層次聚類分析，再與主成分分析得到的結(jié)果綜合分析后，得到了13 個主要變量：氫油比、反應(yīng)過濾器壓差、反應(yīng)器上部溫度、反應(yīng)器底部溫度、烯烴含量、芳烴含量、反應(yīng)器頂?shù)蛪翰?、待生吸附劑性質(zhì)（焦炭含量）、待生吸附劑性質(zhì)（硫含量）、再生吸附劑性質(zhì)（焦炭含量）、再生吸附劑性質(zhì)（硫含量）、辛烷值和硫含量。結(jié)合工藝要求和操作經(jīng)驗發(fā)現(xiàn)這13 個主要變量符合實際情況，并在后面的辛烷值損失模型建立和優(yōu)化過程中，可呈現(xiàn)較好的效果，證明了這13 個主要變量具有有效性和代表性。

2 基于BP 神經(jīng)網(wǎng)絡(luò)的辛烷值損失預(yù)測模型

反向傳播網(wǎng)絡(luò)可以擬合任何非線性函數(shù)，是可以實現(xiàn)從輸入到輸出的端到端網(wǎng)絡(luò)，而且訓(xùn)練網(wǎng)絡(luò)的方法非常簡單、有效，只需要學(xué)習(xí)網(wǎng)絡(luò)參數(shù)即可。通過數(shù)據(jù)預(yù)處理與數(shù)據(jù)降維，將13 個主要變量作為神經(jīng)網(wǎng)絡(luò)的輸入進(jìn)行訓(xùn)練?；贐P 神經(jīng)網(wǎng)絡(luò)的辛烷值損失預(yù)測模型網(wǎng)絡(luò)結(jié)構(gòu)圖如圖1 所示[8]。神經(jīng)網(wǎng)絡(luò)模型可分為兩部分：第一部分神經(jīng)網(wǎng)絡(luò)有800 層，作為特征提取層；另一部分為單層全連接層，作為預(yù)測層。需要注意的是，在模型選擇時，輸入輸出的維數(shù)固定，中間的隱藏層部分需要考慮是否全連接。圖中的箭頭表示數(shù)據(jù)流動的方向，即模型的流程；w表示神經(jīng)網(wǎng)絡(luò)中的權(quán)重矩陣，權(quán)重可通過訓(xùn)練得到；B表示各層神經(jīng)元的偏置，可以使模型具有更魯棒的表達(dá)能力。

圖1 基于BP 神經(jīng)網(wǎng)絡(luò)的辛烷值損失預(yù)測模型網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.1 Network structure diagram of octane loss prediction model based on BP neural network

模型（反向傳播梯度）在更新網(wǎng)絡(luò)參數(shù)時，為避免波動過大出現(xiàn)梯度消失或梯度爆炸現(xiàn)象，導(dǎo)致模型崩潰，需要對輸入進(jìn)行零均值化與歸一化，使輸入數(shù)據(jù)的分布統(tǒng)一。經(jīng)過歸一化，輸入被約束在[-1, 1]的區(qū)間內(nèi)。

第一步零均值化

式中：xi為某一主要變量的某條樣本采樣值。

m為主要變量的樣本個數(shù)。

第二步歸一化方差

在神經(jīng)網(wǎng)絡(luò)中添加激活函數(shù)，可以增加網(wǎng)絡(luò)的非線性表達(dá)能力，并解決線性模型不能很好描述的函數(shù)

關(guān)系[5]。神經(jīng)網(wǎng)絡(luò)中所采用Relu 激活函數(shù)為

式中：z為激活函數(shù)的輸入，在本文的網(wǎng)絡(luò)中為激活函數(shù)前一層神經(jīng)元的輸出；

B為對應(yīng)層神經(jīng)元的偏置（參見圖1）。

Relu 函數(shù)在正區(qū)間解決了梯度消失的問題，而且使神經(jīng)網(wǎng)絡(luò)的收斂速度和計算速度較使用其他激活函數(shù)（如sigmoid 等）有明顯提高，起到了類似Dropout 的各層神經(jīng)元隨機(jī)連接的效果。

考慮到神經(jīng)網(wǎng)絡(luò)的輸入數(shù)據(jù)為篩選出的主要操作變量，故采用均方誤差（mean square error，MSE）來進(jìn)行模型的評價。MSE 可以評價數(shù)據(jù)的變化程度，其值越小說明模型預(yù)測結(jié)果具有更好的精確度[6]，且具有較小的波動。具體表達(dá)式如下：

式中：yi、分別為辛烷值損失真實值和模型預(yù)測值；wi為神經(jīng)網(wǎng)絡(luò)中各預(yù)測值的權(quán)重。

3 模型求解

在對神經(jīng)網(wǎng)絡(luò)模型的實際訓(xùn)練過程中，采用Adam 優(yōu)化器進(jìn)行梯度更新。相比于傳統(tǒng)的隨機(jī)梯度下降，Adam 通過計算梯度矩陣為各個參數(shù)設(shè)置更加靈活的學(xué)習(xí)速率，對變化頻繁的參數(shù)采用大步長進(jìn)行學(xué)習(xí)，對于稀疏的參數(shù)采用小步長進(jìn)行學(xué)習(xí)。第k步梯度gk更新如式（11）所示。

式中α為調(diào)節(jié)系數(shù)，用于調(diào)節(jié)變化快慢。

將整個數(shù)據(jù)集劃分為80%的訓(xùn)練集和20%的測試集。測試集不參與模型訓(xùn)練，在測試集上通過模型的預(yù)測結(jié)果來觀察模型的實際泛化能力。經(jīng)過2 048 Epochs 訓(xùn)練，模型的損失函數(shù)變化如圖2 所示。

由圖2 可知，隨著訓(xùn)練次數(shù)的增加，在測試集上的損失值逐漸降低，在2 000 次左右達(dá)到了9.928×10-6，這說明模型有良好的泛化性能。

圖2 損失函數(shù)的變化圖Fig.2 Change diagram of the loss function

以損失函數(shù)的變化作為本研究中神經(jīng)網(wǎng)絡(luò)預(yù)測模型的測試標(biāo)準(zhǔn)之一。使用Matlab 軟件建立網(wǎng)絡(luò)的訓(xùn)練參數(shù)如下：迭代次數(shù)為2 048，學(xué)習(xí)率為0.01。訓(xùn)練時模型參數(shù)的變化如圖3所示。由圖3a可以看出，隨著epoch 的增加，模型梯度逐漸減小，變化趨勢逐漸穩(wěn)定。圖3b 可以看出在訓(xùn)練時沒有出現(xiàn)模型崩潰現(xiàn)象，具有很好的穩(wěn)定性。

圖3 訓(xùn)練時神經(jīng)網(wǎng)絡(luò)模型參數(shù)變化圖Fig.3 Neural network model parameter variation diagram during training

經(jīng)過訓(xùn)練，訓(xùn)練模型的最終回歸結(jié)果如圖4 所示。由圖4 可以看出，目標(biāo)值和輸出結(jié)果基本上在直線y=x附近分布，說明訓(xùn)練結(jié)果比較理想。

圖4 模型的回歸結(jié)果Fig.4 Regression results of the model

圖5 為各樣本的辛烷值損失實際值與預(yù)測值的分布情況。由圖5 可以看出，模型擬合效果明顯，可以較好地表達(dá)辛烷值的變化情況，并且能取得很好的預(yù)測效果。

圖5 各樣本預(yù)測值與真實值分布圖Fig.5 Predicted value and true value distribution ofeach individual sample

圖6 為3 509 次迭代的損失函數(shù)變化圖。由圖可以看出，在訓(xùn)練2 000 次以后，損失函數(shù)下降趨勢基本趨于穩(wěn)定，與圖2 中2 048 次迭代的結(jié)果相比，MSE在預(yù)測和訓(xùn)練上的變化區(qū)別不大，擬合程度較高，這說明模型具有良好的有效性和泛化能力。

圖6 3 509 次迭代損失函數(shù)變化圖Fig.6 Loss function variation diagram for 3 509 iterations

4 基于遺傳算法的方案優(yōu)化

遺傳算法是一種隨機(jī)化搜索方法，是一種模擬達(dá)爾文進(jìn)化論和孟德爾遺傳學(xué)機(jī)理的計算模型。遺傳算法由編碼、適應(yīng)度評估和遺傳運算3 部分組成[7]。根據(jù)損失散點圖，將損失變化通過三角函數(shù)進(jìn)行擬合，即可得到每一個樣本的損失值。由于每一個樣本包含13 個主要操作變量，損失值可以由操作變量表達(dá)，在此可將問題轉(zhuǎn)變?yōu)榛谥饕僮髯兞康那€優(yōu)化問題。利用遺傳算法解決多元非線性回歸問題，并向最小化損失值的方向優(yōu)化，得到13 個主要操作變量的對應(yīng)取值。根據(jù)變量的取值范圍，得到對應(yīng)主要變量的優(yōu)化操作條件。

對于本問題，使用傳統(tǒng)方法難以獲得最優(yōu)解，而多峰最優(yōu)化問題具有多個局部解的特性，因此使用多峰非線性方法來求解，具體實驗參數(shù)設(shè)置如表4 所示。其中，自變量上下界用來表示每個主要操作變量取值范圍，交叉概率、變異概率用于控制個體改變頻率，設(shè)置迭代次數(shù)為1 000，使算法能夠得到足夠穩(wěn)定的優(yōu)化解。

表4 實驗參數(shù)設(shè)置Table 4 Experimental parameter setting

經(jīng)過實驗，遺傳算法的損失變化如圖7 所示。由圖7 可以看出，隨著迭代次數(shù)增加，損失的波動幅度逐漸減小，最終解逐漸趨于穩(wěn)定。損失變化曲線表達(dá)式為

圖7 遺傳算法的損失變化曲線Fig.7 Loss change curve of genetic algorithm

遺傳算法最優(yōu)個體得分曲線如圖8 所示。由圖8可能看出，遺傳算法在400 次迭代后每代最優(yōu)函數(shù)值變化逐漸穩(wěn)定，這證明算法已逐漸收斂于最優(yōu)解。

圖8 每代最優(yōu)個體得分變化曲線Fig.8 Optimal individual score variation for each generation

通過上述遺傳算法，可以得到一批子代，每個子代對應(yīng)一組變量。根據(jù)辛烷值損失值與13 個變量之間的關(guān)系，建立如下?lián)p失值與主要變量間的表達(dá)式：

由于優(yōu)化問題一般是從凸函數(shù)最小化尋找最優(yōu)解，因此采取最小化f值的方式，獲得x變量的取值，x對應(yīng)13個主要操作變量。經(jīng)過遺傳算法的編碼解碼，從遺傳算法中選擇降幅大于30%的樣本，每個樣本對應(yīng)主要變量優(yōu)化后的一種操作條件。符合條件的優(yōu)化方案，按照降幅比降序排列后，結(jié)果如表5 所示。當(dāng)13 個主要變量分別優(yōu)化至表5 首行對應(yīng)數(shù)值時，所能達(dá)到的降幅比為44.822%，是所有得到的優(yōu)化方案中的最優(yōu)解。

表5 符合條件的優(yōu)化方案Table 5 Satisfactory optimized schemes

當(dāng)某一變量逐步優(yōu)化，其他變量不變時，辛烷值增加量與優(yōu)化步數(shù)之間的變化情況如圖9a和9b所示。由于再生吸附劑性質(zhì)硫優(yōu)化迭代步數(shù)較多，因此通過圖9b 展示。由圖9 可以看出，辛烷值增加量的變化與氫油比、反應(yīng)器底部溫度和反應(yīng)器頂?shù)讐翰钭兓收嚓P(guān)；隨著變量反應(yīng)過濾器壓差值越接近最優(yōu)理想狀態(tài)，對辛烷值的影響變?nèi)?；辛烷值增加量變化與反應(yīng)器上部溫度變化呈先負(fù)相關(guān)后正相關(guān)的關(guān)系，且變化趨勢大致相同，而待生吸附劑中的焦炭含量對辛烷值變化影響相反；辛烷值增加量與烯烴含量變化呈不斷抖動的影響關(guān)系，正負(fù)相關(guān)性交替出現(xiàn)；辛烷值增加量與芳烴含量呈近似指數(shù)變化的正相關(guān)關(guān)系，且越接近最優(yōu)狀態(tài)影響越明顯；再生吸附性質(zhì)（焦炭含量）對辛烷值增加量產(chǎn)生不斷交替的正負(fù)相關(guān)性，但總體呈微弱的正相關(guān)性；辛烷值增加量與再生吸附性質(zhì)（硫含量）呈倒二次函數(shù)變化趨勢；原材料對產(chǎn)品辛烷值增加量的影響呈現(xiàn)出負(fù)相關(guān)的特性，且隨著其值逼近最優(yōu)值，對產(chǎn)品辛烷值增加量的負(fù)向影響力度愈加強(qiáng)烈。由圖9c 可知，當(dāng)所有主要變量同時進(jìn)行優(yōu)化時，主要變量對辛烷值增加量的綜合影響可以分為兩個階段，從一開始的初始狀態(tài)到第10 次迭代，對辛烷值增加量的變化影響明顯；第10 次迭代之后，對辛烷值的影響雖然仍然是正相關(guān)，但其影響程度下降了很多，其中一個原因是有些變量已經(jīng)達(dá)到最優(yōu)解。

圖9 各主要變量變化及辛烷值增量變化軌跡圖Fig.9 Trajectories of major variables with their octane numbers

5 結(jié)語

本文通過層次聚類的方法，對預(yù)處理后的原數(shù)據(jù)集進(jìn)行降維處理，建立基于BP 神經(jīng)網(wǎng)絡(luò)的辛烷值損失預(yù)測模型，對辛烷值損失及其指標(biāo)進(jìn)行預(yù)測；使用遺傳算法迭代得出辛烷值損失降幅比大于30%的優(yōu)化方案及最優(yōu)優(yōu)化策略，并通過可視化將各個操作變量的變化軌跡展示出來加以分析。實驗結(jié)果表明，使用非線性層次聚類法比使用線性聚類法效果更好，能夠很好地對數(shù)據(jù)集進(jìn)行降維，降低了神經(jīng)網(wǎng)絡(luò)模型的訓(xùn)練與收斂難度。預(yù)測模型經(jīng)過2 048 次迭代后，神經(jīng)網(wǎng)絡(luò)模型的損失函數(shù)下降趨勢基本穩(wěn)定。在數(shù)據(jù)集上，采用8:2 的比例構(gòu)建訓(xùn)練集和測試集，并在測試集上達(dá)到了9.982 6×10-6的誤差精度，能滿足模型預(yù)測的要求和精度。根據(jù)樣本硫含量不大于5 μg/g、主要操作變量的取值范圍，以及13 個主要變量的變化步長，利用遺傳算法迭代得到最優(yōu)的參數(shù)集合。從集合中找到子代滿足辛烷值損失降幅比達(dá)到30%的樣本并排序，得到全部優(yōu)化方案以及最優(yōu)優(yōu)化方案。在最優(yōu)方案中辛烷值損失的降幅比達(dá)到44.822%，且各主要變量最優(yōu)值均處于合理取值范圍內(nèi)。結(jié)果表明，本文所采用的數(shù)學(xué)模型與優(yōu)化算法可以很好地對辛烷值損失問題進(jìn)行預(yù)測與優(yōu)化。