基于自適應采樣算法的芳烴異構化代理模型

謝雨珩，李智，楊明磊，杜文莉

（化學工程聯合國家重點實驗室，華東理工大學化工過程先進控制和優化技術教育部重點實驗室，上海200237）

引 言

異構化裝置是芳烴生產過程中的重要環節。隨著芳烴反應動力學模型的深入研究[1]，以及計算機過程模擬技術的發展，研究人員可以利用過程模擬軟件例如Aspen Hysys 對異構化環節構建機理模型，通過改變溫度、壓力、流量等操作參數獲得關鍵產物收率，從而對操作參數進行實時的控制和優化。但是機理模型本身計算時間較長，實際過程中直接使用機理模型會耗費大量計算時間，特別是在進行優化時，由于需要反復調用機理模型，所以往往出現優化時間過長、效率過低等問題，很大程度上限制了機理模型在實際優化中的應用。

代理模型的出現，成為解決這一問題的關鍵技術之一[2]。代理模型是在保證一定精度的情況下，使用計算復雜度較小的數學模型對計算復雜度較大的原始機理模型進行模擬。一方面保證了模型的精度，另一方面大幅提高了計算效率[3]。常用代理模型有：多項式響應面模型[4-5]、徑向基函 數 模 型[6-7]、Kriging 模 型[8-9]、支 持 向 量 回 歸 模型[10]、多元自適應樣條回歸模型[11-12]、人工神經網絡[13-14]等。

代理模型的建立，主要依賴于采樣點的選擇，常用的一次采樣方法包括均勻采樣（uniform design）[15]、拉 丁 超 立 方 采 樣（Latin hypercube sampling，LHS）[16]、增量拉丁超立方采樣（inherited Latin hypercube sampling）[17]等。但是一次采樣方法只是使得樣本點在采樣空間內盡量均勻分布，而實際過程中，模型的某些區域較為復雜，需要較多采樣點，而有些區域實際響應值較為平緩，并不需要過多的采樣點。代理模型與實際模型難免有誤差，必須采取合理的模型管理策略來保證模型的準確性[18]，所以如何使用較少的原始適應度函數評估次數得到機理模型的最佳采樣點分布，進而建立更準確的代理模型成了研究的主要方向。國內外學者提出了很多自適應采樣的方法，Wang 等[19]利用基于委員會的主動學習尋找各子代理模型相差最大的點作為新采樣點，Eason 等[20]提出了基于交叉驗證誤差和最小距離的混合自適應采樣方法，Wei 等[21]提出了基于最大曲率和最小距離的順序采樣方法（MCMPD），Garud 等[22]提出了基于樣本擁擠度和偏離函數的智能采樣算法。Müller 等[23]在當前最優點附近生成新的候選點，對這些候選點的距離準則和響應面準則進行加權獲得新采樣點。與LHS 等一次性采樣方法相比，自適應采樣算法的計算效率更高，采樣點分布更加合理，構建出的代理模型精度也更高。實際應用過程中，由于計算維度往往較高，已有的自適應算法在建模時效率仍然較低，因此需要尋找新的適用于實際大型應用的自適應采樣算法。

為了提高構建代理模型的計算效率，減少原始適應度函數評估次數，本文提出了一種新的計算稀疏度的方法，并結合最鄰近期望提出了一種基于稀疏度和最鄰近期望的自適應代理模型采樣算法（adaptive surrogate sampling algorithm-based sparsity and nearest neighbor，ASSA-SNN）。該算法能夠在模型變化較大的區域以及未探索區域自適應的分配采樣點，從而使得代理模型的精度得到顯著的提高。使用多個測試函數驗證了該方法的準確性和高效性，并將其應用于芳烴異構化環節建模中，實現對異構化產物的準確預測。

1 Kriging代理模型

Kriging 模型是對空間分布的數據求線性最優、無偏內插估計的一種插值方法，被廣泛用于非線性問題代理模型的構造。最早于1951 年由南非地質學家Krige[8]提出，早期主要用來預測礦產儲量的分布，之后Sacks 等[24]將Kriging 理論進一步推廣到確定性計算機實驗的設計和分析領域，并給出了一種較實用的Kriging 模型。Kriging 模型近些年來在化工領域也有較多應用，例如苯乙烯流程優化[25]以及碳二加氫反應器建模[26]。

已知n維的訓練樣本X=[x1,x2,…,xN]T對應的響應值為Y=[y1,y2,…,yN]T。那么，Kriging 的y(x)由一個全局模型和一個局部偏差構成

式中，f(x)是確定性部分；Z(x)是隨機過程，并滿足統計特性：期望為0，方差為σ2，協方差為σ2R(xi,xj)。其中，R為相關函數，通常采用Gauss函數

其中，n為X的維數，θk是待定的模型參數。任意x處的Kriging 模型預測值y?(x)和預測方差s?2表示為

其中，F=[1,1,…,1]T，rT(x)為待預測點x和訓練樣本點X=[x1,x2,…,xN]T之間的相關向量，定義如下

2 基于稀疏度與最近鄰期望的自適應采樣算法

2.1 稀疏度

代理模型的建立需要在機理模型上獲得初始樣本集X，然后根據已有信息找到下一個采樣點。新采樣點需要同時滿足兩個條件：能反映復雜區域的響應值和未搜索區域的響應值。所以需要綜合考慮新采樣點的位置信息和模型響應值信息，前者主要負責全局搜索，后者負責局部搜索。

對于新采樣點的位置信息，本文提出了稀疏度的概念。根據新采樣點xnew與已有采樣點的歐式距離以及每一維數據的大小排序，得到在新采樣點在整個搜索空間的稀疏度。已有采樣點X={x1,x2, …,xN}，搜索空間的上下限分別為LU={lu1,lu2,…,luD}和LD={ld1,ld2,…,ldD}，D為樣本維度。搜索空間內的任意一點xnew的稀疏度Sparsity(xnew)定義如式(7)所示，具體流程如圖1所示。

當函數為2維時，樣本點的稀疏度是一個面積，以圖2 為例，在搜索空間x1,x2∈[-4,4]區間內先隨機生成20 個采樣點。矩形表示[-2.5,-2.5]、[-2.5,3]和[3,3]三個點的稀疏度。明顯看出，[-2.5,3]附近采樣點較多，所以稀疏度較小，而[-2.5, -2.5]和[2.5,3]附近點較少，稀疏度比較大。

當輸入樣本為三維時，稀疏度就可以類比成三維體積，輸入樣本為四維時稀疏度為四維體積。傳統的使用歐式距離計算新采樣點與已有采樣點的位置信息時，往往只考慮了新采樣點與最近已有采樣點的位置信息，忽略了它與其他相鄰采樣點的信息。本文中的稀疏度則綜合考慮了新采樣點附近已有采樣點的位置信息。稀疏度越大，說明這個點附近的采樣點較少，應該著重搜索。稀疏度越小，說明這個點附近采樣點已經很多，不再需要過多的采樣點了。但是，如果只考慮稀疏度就忽略了局部復雜區域的信息，例如最優值附近往往梯度較大，應該需要更多的采樣點。所以這里還需要考慮代理模型的響應值。

2.2 最鄰近期望

代理模型的建立，不僅需要考慮全局信息，還需要對模型的關鍵區域進行重點采樣，模型的關鍵區域往往梯度較大，本文采用最鄰近期望（nearest neighbor expected,NNE）來近似梯度。

已有采樣點X={x1,x2,…,xN}，對應的實際輸出為Y={y1,y2,…,yN}，根據X和Y構建的代理模型為y?。最鄰近期望NNE為

ynearest為樣本集X={x1,x2,…,xN}中距離xnew最近的樣本點的對應的真實響應值。NNE 越大說明新采樣點附近的近似梯度較大，即函數波動較大，需要著重進行采樣。傳統的交叉驗證的主要優點是不需要新的樣本點來衡量構建的代理模型的準確性[27]，但是會消耗大量計算時間，而NNE既可以反映函數的局部信息，還可以大量縮短計算時間，尤其是在維數較高的情況下，交叉驗證往往需要耗費大量的計算時間，而這時，使用最鄰近期望進行局部搜索可以減少大量不必要的計算，使得搜索更有效率。

2.3 自適應采樣策略

圖1 稀疏度計算流程圖Fig.1 Flow chart of sparsity calculation

代理模型建模過程中，最主要的問題是如何利用盡量少的原始函數評估次數，使得樣本點盡可能多地分布在具有關鍵信息的區域[28]，從而縮短總體建模時間。即樣本點個數相同的情況下，使獲取的關于未知設計空間的信息最大化[29]，提升代理模型精度。

利用自適應采樣策略構建代理模型的通用過程步驟如下：

（1）選用適當的取樣方法生成初始樣本點集，并求得其對應的原始適應度函數。

（2）基于上述初始樣本點及其原始適應度函數，建立目標問題的代理模型。

（3）根據已有樣本和已構建的代理模型，尋找能反映原始函數關鍵信息的新采樣點，評估其原始適應度函數并添加到樣本集中。

（4）判斷是否符合終止條件，符合終止條件則輸出已有的代理模型，否則返回步驟（2）。

整個流程的關鍵就是步驟（3）中如何找到反映原始函數關鍵信息的采樣點。在搜索新采樣點過程中代理模型應與原始適應度函數一起使用[30]。新采樣點一方面要考慮全局信息，即新采樣點應該在采樣點稀疏的空間內，以保證代理模型在整個搜索空間的泛化性，同時要考慮響應值的局部信息，即波動較大的區域應該有較多的采樣點以保證代理模型在此區域的準確性，本文根據稀疏度和最鄰近期望的乘積最大來尋找新采樣點。

圖2 二維稀疏度示意圖Fig.2 Diagram of 2D sparsity

稀疏度Sparsity負責控制全局搜索，最鄰近期望NNE負責局部關鍵信息搜索。R為采樣空間定義域。初始樣本點通過拉丁超立方采樣獲得，終止條件為新增采樣點個數達到設定值。使用模式搜索法[31]求解式(9)?？紤]到搜索過程中需要計算歐式距離，但是實際應用中，很多情況下不同維的數據的量綱不同時，歐式距離可能會有較大誤差，所以在構建代理模型前，將所有數據根據采樣空間的上下限進行歸一化，這樣保證歐式距離的計算不會受到量綱的影響。具體流程如圖3所示。

基于稀疏度與最近鄰期望的自適應采樣算法（ASSA-SNN）在尋找新采樣點時可以有效地權衡全局搜索和局部搜索，從而在迭代過程中尋找到對代理模型精度影響較大的點，不斷充實樣本集，提高代理模型的精度。

3 測試函數驗證

本文首先以二維PK 函數為例，PK 函數表達式為

圖3 算法流程圖Fig.3 Flowchart of algorithm

初始采樣點20，新增采樣點80，通過ASSASNN 得到的樣本點分布如圖4 PK 函數采樣點分布情況所示，可以看出在新增的采樣點主要集中在上下兩個極值點附近，而其余較為平坦的區域采樣點較少。使用均方根誤差（root mean square error，RMSE）來對算法的精度進行評估，RMSE表達式為

圖4 PK函數采樣點分布情況Fig.4 Sampling point distribution of PK function

圖5 PK函數RMSE變化情況Fig.5 Variation curve of RMSE of PK function

表1 PK函數自適應采樣與拉丁超立方采樣結果對比Table 1 Comparison of PK function adaptive sampling and Latin hypercube sampling results

本文共選取了8 個測試函數，并與混合自適應采樣算法[20]（MASA）、基于最近鄰和馬氏距離的自適應采樣算法[32]（ASA-NNMD）、智能采樣算法[22]（SSA）進行對比。為了保證對比的有效性，統一使用Kriging 模型，其中，f(x)為零階多項式回歸函數，相關函數為高斯函數[式(2)]。初始采樣點均使用拉丁超立方（LHS）采樣獲得。由于模型的維數越高意味著模型的復雜程度增加，所以針對不同維數的模型，所選取的初始采樣點個數與新增采樣點個數不同，具體采樣方式如表2所示。

為消除不同初始采樣點以及算法尋優時的隨機性，每種算法在每個測試函數上運行10 次取平均值以消除隨機性帶來的影響。最終結果如表3 所示。從結果中可以看出，8 個測試函數中，ASSA-SNN 有6 個測試函數的RMSE 占優，Hartman3 函數的結果也與最小值相差較小，只有在AlpineN.2 函數的結果不如MASA 和SSA。這是由于該函數的極值點過多，局部波動較大，導致搜索到的點過于集中在各個極值附近，影響了整體的效率。

表2 采樣方式Table 2 Sampling configurations

表3 各測試函數的平均RMSETable 3 Average RMSE obtained by four sampling approaches for test cases

為了更好地說明算法的計算效率，表4 列出了每種測試函數在不同方法下運行的平均時間。從表格中可以看出由于SSA 使用了逐一交叉驗證，MASA 使用了K 折交叉驗證，導致計算時間隨著維數增加而迅速上升。ASA-NNMD 由于需要通過計算大量矩陣計算馬氏距離，以及本身的優化過程較為耗時，導致計算時間也較長。而ASSA-SNN 只有基本的歐式距離計算，通過最鄰近期望代替交叉驗證，使用模式搜索代替耗時的智能優化算法，達到了精度與效率的雙贏，不僅精度有提高，計算時間也大幅縮短。

綜上所述，本文提出的ASSA-SNN 在大部分測試函數中的表現優于其他兩種算法，并且維數大于3 后ASSA-SNN 計算時間要明顯優于其他的算法，函數測試結果說明了ASSA-SNN 的建模精度和建模效率要優于其他算法。

圖6 異構化過程流程圖Fig.6 Flow chart of isomerization process

表4 各測試函數的平均計算時間/sTable 4 Average running time obtained by four sampling approaches for test cases/s

4 芳烴異構化代理模型構建實例

芳烴異構化反應的目的是將異構體鄰二甲苯(OX)、間二甲苯(MX)和乙苯(EB)轉化為價值更高的對二甲苯(PX)[33]。異構化過程主要由反應器和分離系統組成，其工藝流程如圖6 所示，首先是反應器，其主要功能就是在催化劑作用下進行異構化反應，將含有貧PX 的C8芳烴混合物轉化成接近熱力學平衡的C8芳烴混合物。其次是分離系統，反應產物在進行分離前需要實現氣液分離，將循環氫和反應液送入不同的裝置。分離的液體經換熱后送入脫輕組分精餾塔，實現輕組分和重組分的分離，分離出來的C8及C8+送到二甲苯精餾單元的分離塔實現C8芳烴與重芳烴分離。

采用化工過程動態模擬軟件Aspen Hysys 建立其機理模型。由于機理模型的計算采用序貫模塊法逐個裝置求解，優化時每次迭代都要將優化后的輸入數據重新送入模型中并調用模型求解，采用種群數為20 的遺傳算法，迭代次數80 次，總共耗時13059.68 s，耗時較長，不利于優化進行，所以尋求利用代理模型節約計算時間。根據異構化環節的反應機理和經驗知識，本文選取8個輸入變量：異構化進料、循環氫流量、補充氫流量、異構化反應溫度、異構化反應壓力、乙苯含量、MX 含量和OX 含量。篇幅所限，本文只選取兩個重要的質量指標進行說明：輕烴收率和C8+收率。利用拉丁超立方采樣獲得80個初始采樣點和1000個測試樣本點，算法最大迭代次數為320，分別用四種方法建立代理模型并計算它們的RMSE。結果如表5所示，所用建模時間如表6所示。

表5 芳烴異構化重要指標的RMSETable 5 RMSE of important indicators for isomerization of aromatics

圖7 給出了異構化代理模型中MX 含量、OX 含量對出口C8+收率的影響關系，可以明顯地反映出模型變量之間的復雜映射關系。圖8 給出了MX含量、OX 含量樣本分布情況，可以看出在各極值附近均有較多樣本點分布。從表5 和表6 看出，四種關鍵性能指標的建模中，ASSA-SNN 的準確度高于其他三種算法，輕烴收率的RMSE 分別降低了16.49%、5.56%、8.20%，輕烴收率的計算時間分別 縮 短 了33.92%、52.41%、89.27%，C8+收 率 的RMSE 分別降低了18.91%、12.20%、7.07%，C8+收率的計算時間分別縮短了35.94%、53.60%、89.61%。以輕烴和C8+收率RMSE 變化曲線如圖9 和圖10 所示，可以看出ASSA-SNN 的整體RMSE 下降速率要高于其他三種方法，最終的RMSE 也小于其他方法。

表6 芳烴異構化重要指標的建模時間/sTable 6 Modeling time for important indicators of aromatic isomerization/s

圖7 異構化代理模型Fig.7 Surrogate model of isomerization

圖8 MX含量、OX含量樣本分布Fig.8 Sample distribution of MX content and OX content

圖9 異構化的輕烴收率RMSE變化曲線Fig.9 Variation curves of RMSE of isomerized light hydrocarbon

圖10 異構化C8+收率RMSE變化曲線Fig.10 Variation curves of RMSE of C8+yield with number of iterations

從測試結果中看出，ASSA-SNN 相比于其他三種方法，精度有所提升，能夠更加準確地對各關鍵性能指標進行建模，并且計算時間大幅縮短，這使得整體的建模效率大為提升，為后續的操作變量和操作條件的優化提供便利。

5 結 論

本文提出了一種新的基于Kriging 代理模型的自適應采樣算法，該算法通過最大化稀疏度與最鄰近期望的乘積來確定新采樣點的位置，從而使得模型能夠根據實際函數自適應的分配采樣點位置。相比其他自適應采樣算法，使用稀疏度可以更好地判斷新采樣點與已有采樣點的位置關系，從而更好地進行全局搜索，用最鄰近期望代替交叉驗證，提高了局部搜索效率，使用模式搜索解決最優問題，極大地節約了計算時間，提升了尋優效率。將該算法與其他的自適應采樣算法在8個測試函數上進行了驗證，表明在相同采樣點個數的情況下，ASSASNN 得到的采樣點分布能夠有效提高建立的代理模型精度并且縮短建模時間。芳烴異構化實例也證明了算法的有效性，為后續的異構化環節操作條件的優化提供了便利。