基于回聲狀態網絡的煉油廠氫氣網絡動態用氫成本預測及應用

左春玲， 張正明， 曹萃文

(華東理工大學 化工過程先進控制和優化技術教育部重點實驗室，上海 200237)

近年來，隨著煉油產品的升級和原油的重質化、劣質化，石油化工企業對氫氣資源的需求迅速增加，煉油企業催化重整裝置的副產氫氣已經滿足不了生產需要，需采取增加制氫裝置、擴大已有裝置制氫能力、從外界購買新氫、優化氫氣網絡等多項措施來彌補不足。氫氣從產出到被下游裝置使用，如何合理地安排好產、耗平衡，直接涉及到煉油廠整體的平穩生產和經濟效益。

已有的對氫氣優化利用的研究，主要為夾點分析法[1-3]和超結構法[4-9]。夾點分析法是對現有裝置的氫氣純度、壓力、流量、雜質含量等約束條件進行分析[10]，主要以公用工程最小用量[2]或最小新氫供應量[3]為目標，通過繪制氫負荷流量圖等來確定夾點位置。該方法易于實際工程操作和執行，方便確定網絡整體用氫目標，但不適合進行整個系統的動態經濟效益指標分析。

超結構法[4]簡化了裝置，注重各裝置之間的關系，將產氫、耗氫裝置分別作為氫源、氫阱，建立了包括所有可能有連接關系的網絡超結構。這種連接關系考慮實際生產過程中的各種約束條件，如：物料守恒、濃度守恒、壓力條件和設備配置等，用數學規劃法求解，滿足所有約束條件并使目標函數最優。在優化煉油廠的氫氣網絡時，通常以靜態成本最小為目標[8,10]，可以綜合考慮操作費用和設備投資費用。2012年，Liao等[8]提出了一種適用于氫氣網絡改造設計的方法，通過建立狀態空間超結構模型獲得了更多的網絡結構，權衡了各種靜態操作成本和投資成本。曹萃文等[9]發明了一種在超結構法中混合夾點約束的煉油廠氫氣網絡優化調度方法，其中的成本目標仍然是靜態成本。康永波[11]提出以供氫成本和新制氫氣剩余量為目標的氫氣網絡多目標模型，在氫耗需求不確定環境下，建立魯棒優化模型，利用遺傳算法求解，并與普通優化比較，證明了在氫氣網絡優化過程中充分考慮不確定參數變化的重要性。

氫氣網絡是復雜的非線性系統，對復雜系統進行研究主要是建立系統的差分方程模型或者微分方程模型，根據數學模型的特點來探索系統的動力學機制。然而，在實際問題中，由于系統結構復雜，參數隨時間變化，通常無法獲得完備的和確定性的信息，很難建立精確的解析模型。在此情況下，可以運用數據驅動技術研究系統的動態行為。該技術對系統內在的背景知識需求較少，對模型機理要求較少，為那些無法建立精確解析模型的實際系統提供了一條新的研究途徑。2004年，Jaeger提出了回聲狀態網絡(Echo state network，ESN)，并對混沌時間序列進行預測，通過實驗證明回聲狀態網絡的預測精度比之前報道的預測結果提高了2400倍[12]。回聲狀態網絡因為其訓練簡單和建模精度高等優點，廣泛應用于時間序列預測[12-13]、非線性系統識別[14]、圖像處理[15]、自然語言識別[16]、電力負荷預測[17]等領域，針對不同應用領域對回聲狀態網絡的優化改進算法研究也越來越多。田中大等[18]將一種基于貝葉斯回歸的回聲狀態網絡應用于小型水電站的短期電力預測，克服了傳統神經網絡中存在的病態解降低模型泛化能力的缺陷。Deihimi等[19]提出一種多層的回聲狀態網絡結構，用于預測特定時刻的電力負荷。模型可預測一天24 h每個整點時刻的電力負荷，使預測結果更加精準且符合實際應用的要求。王莉莉[20]針對污水處理過程的非線性、大時變等特點，提出一種基于回聲狀態網絡的多變量自適應預測控制策略，實現了優化控制。劉穎等[21]通過奇異值分解求取ESN儲備池到輸出層的權值對ESN進行改進，改進后的ESN已應用于高爐煤氣發生量的預測中。盛春陽等[22]提出一種基于濾波過程改進回聲狀態網絡的時間序列預測方法，采用濾波過程來估計回聲狀態網絡儲備池神經元的狀態和網絡輸出權值，能避免參數求解過程中常出現的異常解問題，有效預測冶金企業氧氣系統的氧氣流量。

筆者以某煉油廠實際焦化干氣水蒸氣制氫裝置產氫量的動態變化數據為基礎，通過“基于LabVIEW的氫氣網絡優化調度軟件V 1.0”[23]計算得到的動態成本數據，詳細分析了某煉油廠氫氣網絡成本特點，首次提出采用預測精度較高的回聲狀態網絡對氫氣網絡動態成本進行數據驅動建模及預測，然后根據預測出的成本在不同工況下最低供氫成本的優化操作區域中的位置，對氫氣網絡的合理配置提供決策支持。

1 回聲狀態網絡

1.1 回聲狀態網絡模型

經典的回聲狀態網絡是一種三層遞歸神經結構[13,19]，包括K個神經元組成的輸入層、N個神經元組成的隱含層(儲備池)和L個神經元組成的輸出層，如圖1所示。隱含層又稱儲備池，含有成百上千個稀疏遞歸連接的神經元，與遞歸神經網絡中的隱含層對應。

ESN狀態方程：

x(t)=f(Winu(t)+Wx(t-1)+Wbacky(t-1))

(1)

ESN輸出方程：

y(t)=fout(WoutX(t))

(2)

公式(1)和(2)中：x(t)∈RN×1、x(t-1)∈RN×1(R為實數集合)分別為t時刻、t-1時刻儲備池內部狀態；u(t)∈RK×1是離散時刻t時的外部輸入；Win∈RN×K為輸入連接權值矩陣；W∈RN×N為儲備池內部神經元連接權值矩陣(一般保持[1%,5%]的稀疏連接)，Win和W均隨機初始化并且保持固定不變；Wback∈RN×L為輸出反饋矩陣(目的是將上一時刻網絡輸出反饋到當前時刻狀態中；X(t)≡[x(t)T,u(t)T]T；Wout∈RL×(N+K)是輸出連接權值矩陣；y(t)∈RL×1、y(t-1)∈RL×1分別是離散時刻t、t-1時對應的ESN輸出。f(·)表示儲備池神經元激活函數，通常取非線性函數，如雙曲正切函數或sigmoid函數；fout(·)表示輸出層神經元激活函數，通常取線性函數，對于高度非線性的對象，也可選用非線性函數。

圖1 標準ESN的結構示意圖Fig.1 Structure of standard ESN

儲備池是回聲狀態網絡的核心結構，其參數對回聲狀態網絡的性能好壞有很大的影響。根據不同時間序列的特點，設計適合的儲備池結構是回聲狀態網絡建模的首要問題。儲備池主要參數[15,18]包括激活函數類型、儲備池的規模、內部連接矩陣的譜半徑、稀疏度以及輸入變換系數。

1.2 回聲狀態網絡的訓練與測試

首先根據實際樣本數據來確定回聲狀態網絡的規模和參數，再進行網絡的訓練和測試。回聲狀態網絡的訓練過程就是根據給定的訓練樣本u(t)、y(t)(t=1,2,…,r)確定系統中的輸出連接權矩陣Wout的過程。一般假定Wback=0。ESN的訓練過程可以分為2個階段：采樣階段和權值計算階段。

1.2.1 采樣階段

采樣階段是首先任意選定網絡初始狀態，但是通常情況下令x(0)=0。訓練樣本u(t)(t=1,2,…,m,…,r)經過輸入連接權矩陣Win，y(t)經過反饋連接權矩陣Wback分別被加入到儲備池。按照公式(1)狀態更新方程，依次完成系統狀態的計算。

為了計算輸出連接權值矩陣Wout，需要從某一時刻m開始收集儲備池狀態變量x(t)和輸入變量u(t)，并以向量[xT(t),uT(t)]=[x1(t),x2(t),…,xN(t),u1(t),u2(t),…,uK(t)](t=m,m+1,…,r)為行構成矩陣H(M-m+1,N+K)，其中M是訓練樣本的數量；同時相應的目標輸出y(t)也被收集，以yT(t)=[y1(t),y2(t),…,yL(t)](t=m,m+1,…,r)為行構成矩陣Z(M-m+1,L)。

1.2.2 權值計算

(3)

(4)

Wout=((HTH)-1HTZ)T

(5)

2 基于回聲狀態網絡的氫氣網絡動態用氫成本預測

煉油企業各類產品的生產均隨市場需求和煉油廠實時工況的變化而變化，因此，對氫氣網絡因產氫量和耗氫量變化以及不同工況導致的氫源(產氫裝置)、氫阱(耗氫裝置)不同配置下的實時動態用氫成本進行準確的預測，能夠對氫氣網絡的合理配置與優化提供決策支持。ESN的儲備池神經元連接方式復雜，能夠與其他神經元隨機連接，也可以自連接，相比于BP、KNN等及其他遞歸神經網絡，可以實現復雜系統的高精度建模。

2.1 數據來源及其處理

筆者對國內某大型煉油廠生產過程進行了調研，該煉油廠常規使用的供氫裝置有外購純氫(氫源1)、2套水蒸汽制氫裝置(氫源2)、2套乙烯裝置附產氫(氫源3)、3套重整裝置附產氫(氫源4)、2套裂解汽油加氫裝置附產氫(氫源5)、1套柴油加氫裝置附產氫(氫源6)、1套催化汽油脫硫裝置附產氫(氫源7)、1套高壓加氫裂化裝置附產氫(氫源8)、2套歧化裝置附產氫(氫源9)和2套異構化裝置附產氫(氫源10)，其某月的氫源數據平均值見表1[11]。該煉油廠耗氫裝置有3套裂解汽油加氫裝置(氫阱1)，3套柴油加氫裝置(氫阱2)，1套催化汽油脫硫裝置(氫阱3)，1套中壓加氫裂化裝置(氫阱4)，1套高壓加氫裂化裝置(氫阱5)，1套航煤臨氫脫硫醇裝置和石腦油預加氫裝置(氫阱6)，1個滌綸部用氫(氫阱7)，1套渣油加氫裝置、塑料部用氫和硫磺裝置(氫阱8)，2套歧化裝置(氫阱9)，2套異構化裝置(氫阱10)，其同月的氫阱數據平均值見表2[11]。(表1和表2中的數據，是排除地理位置因素，將相同氫氣純度的裝置歸為了同一個等級后的結果。)

表1 某煉油廠的氫源數據Table 1 Hydrogen sources data of the refinery

1) Normal state

筆者在“基于LabVIEW的氫氣網絡優化調度軟件V 1.0”[23]中運用了文獻[11]的超結構線性規劃數學模型，列出如下：

(1)模型約束

式(6)為氫源流量約束，表明氫源j分配到所有氫阱k的氫氣流量不大于氫源j的供氣流量，且供應流量必須不大于氫源j的最大設計供應量。

(6)

表2 某煉油廠的氫阱數據Table 2 Hydrogen sinks data of the refinery

1) Normal state

(7)

式(8)為氫氣濃度約束。該式表明各氫源j流入氫阱k的氫氣純度必須不小于氫阱k的氫氣純度要求。

(8)

其中，ωj為氫源j的氫氣純度，%；ωk為氫阱k要求的氫氣純度，%。

(2)目標函數：

超結構模型的目標函數可以設為多種形式，筆者為考察煉油廠氫氣網絡的動態用氫成本，列目標函數如式(9)所示。

(9)

其中，TH2為單位時間H2網絡的耗氫總成本，CNY/h；cj為氫源j的H2價格，CNY/m3。

常規研究一般只針對表1和表2中的H2網絡各裝置產氫和耗氫數據，代入超結構模型進行求解后得到一個固定的網絡用氫成本或其他目標，無法反映實際生產的動態環境對目標的影響。為改善這種情況，筆者以該煉油廠焦化干氣水蒸氣制氫裝置(氫源2)為例，采集了裝置的實際產氫數據，并在Aspen HYSYS軟件中模擬該裝置的生產過程，在軟件的各項輸出數據與實際數據吻合后，擴展主要操作參數變動值為各類不同工況下的極限值，得到4100組氫源2的產氫數據。圖2為對應的不同壓縮機出口壓力和原料氣中水/碳摩爾比條件下的氫源2的產氫量。在筆者開發的應用軟件[23]中，對每一組數據采用單純型算法求解氫氣網絡線性規劃模型得到全局最優解，即得到在氫源2產氫量動態變化下的4100組成本和相應的氫氣網絡配置結構。最后將所有影響產氫量的參數配置數據、氫氣流量和純度數據、成本優化結果都通過LabVIEW軟件自動存入MySQL數據庫中。圖3為通過不同壓縮機出口壓力和原料氣中水/碳摩爾比條件下的最小用氫成本。

圖2 不同壓縮機出口壓力和原料氣水/碳摩爾比條件下的氫源2產氫量Fig.2 Hydrogen production of hydrogen source 2 underdifferent compressor outlet pressure andfeedstock water/carbon molar ratios

2.2 氫氣網絡動態成本預測

由2.1節的數據處理過程可以看到，隨著氫源2的產氫量變化(實際上所有氫源、氫阱的產氫量和耗氫量及相關指標均可能發生變化)，氫氣網絡的實時動態用氫成本并不容易實時在線得到，因此避開氫氣網絡中各裝置機理模型在Aspen HYSYS軟件中的復雜仿真過程，開展根據各臨氫裝置主要操作參數變化導致的氫氣網絡成本變化的數據驅動建模及實時預測，對提高煉油廠的操作優化水平、降低氫氣使用成本具有很好的指導意義。

在焦化干氣水蒸氣制氫裝置(氫源2)的生產過程中，原料氣的流量、水/碳摩爾比、烯烴含量、硫含量、催化劑活性、反應器溫度和壓力的改變都會影響氫氣產量，但在實際生產中，影響產氫量最大的常規操作參數為原料氣的水/碳摩爾比和壓縮機出口壓力數據(壓縮機位于第一級主反應器入口處)。

圖3 最小用氫成本與水/碳摩爾比和壓縮機出口壓力的關系Fig.3 Relationship between minimum hydrogen cost,water/carbon molar ratio and compressor outlet pressure(b) is vertical view of (a)

下面以4100組氫氣網絡的動態用氫成本數據為例，說明訓練回聲狀態網絡，從而預測氫氣網絡的最新動態用氫成本的過程。

焦化干氣水蒸氣制氫裝置的原料氣水/碳摩爾比和壓縮機出口壓力為回聲狀態網絡的二維輸入，所對應的最小用氫成本作為目標輸出。4100組數據分為兩部分，3000組為網絡訓練數據集，1100組為網絡測試數據集。回聲狀態網絡的性能通過使用公式(10)的歸一化均方根誤差(Normalized root mean square error，e)來評估[12]。

(10)

儲備池內部取1500個神經元(根據多次改變儲備池內部神經元個數的實驗結果比較，N=1500能有更好的預測效果)，稀疏連接為2%，隨機生成輸入連接矩陣Win和儲備池內部連接矩陣W并保持不變，譜半徑是儲備池內部連接權矩陣W的最大特征值的絕對值。運用訓練數據集和任意網絡的初始狀態(令x(0)=0)進行儲備池狀態的更新計算，最后求解網絡輸出值，可以得到如圖4、圖5、圖6所示的結果。因儲備池內部的神經元個數較多，圖4顯示的是儲備池內部任意選取的4個神經元的前200個狀態變量值(x1(t),x2(t),x3(t),x4(t),t=1,2,…,200)，被收集在矩陣H中；不同的輸入變量在狀態更新方程(見公式(1))計算后產生不同的狀態變量x(t)。圖5為在某一次獨立運行時3000組訓練數據目標期望值與回聲狀態網絡預測輸出值的對比圖。圖6為1100組測試數據目標期望值與回聲狀態網絡預測輸出值的對比圖。由圖5和圖6可以看出，在假設其他條件都不變的情況下，回聲狀態網絡可以很好地根據壓縮機出口壓力和原料氣水/碳摩爾比來預測氫氣網絡的最小用氫成本。

圖4 選取的4個儲備池內部神經元的前200個狀態變量值Fig.4 The front 200 state variable values for neurons in the four selected reserve pools(a) X1(t); (b) X2(t); (c) X3(t); (d) X4(t)

表3為回聲狀態網絡在氫氣網絡動態成本預測中的性能指標。性能評估指標為歸一化均方根誤差、平均絕對誤差、相對誤差小于1%、1%～5%和大于5%的樣本比例，結果是獨立運行50次計算得到的平均值。由表3可以看出，多于95%的數據可以達到預測相對誤差小于5%。這足以區分該條件下的最小用氫成本所在的范圍，并為最終能尋找到最優氫氣網絡結構和操作條件提供數據支撐。

3 基于氫氣網絡動態用氫成本預測值的區塊操作優化

3.1 區塊劃分分析

為將第2.2節中預測出的實時網絡用氫成本數據用于指導運行操作，筆者對氫源2的2個典型操作參數(原料氣的水/碳摩爾比和壓縮機出口壓力)變化時產生的4100組用氫成本數據進行了分析。由圖3 可知，最小用氫成本與2個典型參數之間的關系是非線性的，可以以壓力2840 kPa為分界點將其劃分為12個區塊：圖7(a)顯示了壓力小于等于2840 kPa時最小供氫裝置成本與水/碳摩爾比和壓縮機出口壓力的關系(6個區塊)；圖7(b)顯示了壓力大于等于2840 kPa時最小供氫裝置成本與水/碳摩爾比和壓縮機出口壓力的關系(6個區塊)。

3.2 在區塊跳躍處的操作優化分析

當最小用氫成本出現區塊間跳躍時(如圖7所示不同顏色的區塊)，說明該最小用氫成本對應的氫氣網絡中氫源、氫阱的氫氣流量分配發生了改變。隨機選取其中一個跳躍的相鄰部分，如圖7(a)圓圈圈出來紅色區域和綠色區域的兩點。紅色區域上的圓圈(壓力為2790 kPa、水/碳摩爾比為3.126)，此時各氫源分配到各氫阱之間的氫氣流量如表4所示，最小用氫成本為287851.70 CNY/h；綠色區域上的圓圈(壓力為2790 kPa、水/碳摩爾比為3.165)，此時各氫源分配到各氫阱之間的氫氣流量如表5所示，此時的最小用氫成本為255635.49 CNY/h。從表4和表5可以看出，由于水/碳摩爾比(即操作方式)的微小變化，導致氫氣網絡中氫源4、氫源5和氫源8對氫阱1、氫阱2、氫阱3和氫阱4的供氫網絡結構和氫氣流量發生改變，以至于最小用氫成本減少32216.21 CNY/h。因此，當第2.2節中預測的最新網絡用氫成本數據處在區塊的邊界處時，對應的操作參數應向成本較小的下區塊對應值處調整。

圖5 訓練數據的期望輸出值和ESN輸出值Fig.5 Expected value and ESN output value of training data(a) All of the training data； (b) Part of the training data Teacher sequence; Predicted sequence

圖6 測試數據集的期望值和ESN輸出值Fig.6 Expected values and ESN output values of test data(a) All of the test data； (b) Part of the test data Teacher sequence; Predicted sequence

ItemeAverageabsoluteerrorRatio of relative error/%<1%1%-5%Train0.3013321.44867.926.7Prediction0.3123431.08267.126.9

3.3 在同一區塊上的操作優化分析

由圖7可知，在一定的水/碳摩爾比和壓縮機出口壓力范圍內(圖7(a)和圖7(b)中相同顏色區域，即同一區塊上)，氫氣網絡結構不會發生頻繁的變動，體現了實際可操作性。即：如果第2.2節中預測的網絡最新用氫成本值在某一區塊中時，可以根據實際情況調節水/碳摩爾比或壓力到同一區塊的最小用氫成本時的條件。圖7中三角形標出的點即為每個不同區塊的最小用氫成本值。如:當實際需要水/碳摩爾比范圍在3.946～5.079、壓力范圍在2840～3270 kPa時(圖7(b)中最下方深藍色區域)，操作調節水/碳摩爾比為4.063和壓力為3010 kPa，則可以在不改變氫氣網絡結構的前提下，獲得最小用氫成本為 248250.00 CNY/h。

表4 氫源-氫阱之間的氫氣流量分配(圖7(a)紅色區域的圓圈)Table 4 Flow distribution of hydrogen between hydrogen sources and sinks (Circle in the red region of Fig.7(a)) QH2/(m3·h-1)

The flow rate distributions of Hydrogen source 3, Hydrogen source 7, Hydrogen source 9, Hydrogen source 10 and Hydrogen sink7, Hydrogen sink 8, Hydrogen sink 9, Hydrogen sink10 are all 0, so they are not listed in the table.

圖7 最小用氫成本與水/碳摩爾比和壓縮機出口壓力的關系Fig.7 Relationship between the minimum hydrogen cost, the water/carbon molar ratio and compressor outlet pressure(a) p≤2840 kPa; (b) p≥2840 kPaTriangles mark the state, in which each block costs the least. Circles are examples of optimization operations.

QH2/(m3·h-1)

Same legend as in Table 4.

4 結 論

(1)在利用某煉油廠氫氣網絡實際數據和Aspen HYSYS流程模擬獲得的數據集合基礎上，分析了氫氣網絡最小用氫成本和相應的氫氣網絡流量分配變化的情況。

(2)以該煉油廠氫氣網絡中水蒸汽制氫裝置的原料氣水/碳摩爾比和壓縮機出口壓力數據變化對最小用氫成本的影響為例，說明了運用回聲狀態網絡對氫氣網絡最小用氫成本進行預測的可操作性，超過95%預測結果的相對誤差小于5%。

(3)由于最低用氫成本與2個主要操作參數(原料氣的水/碳摩爾比和壓縮機出口壓力)之間存在非線性關系，對其進行了區塊操作優化劃分分析。結果表明：可以根據預測的最小用氫成本在不同區塊上的位置，在對應實際生產條件范圍內尋找到最低用氫成本時的最優操作參數設置；也可以根據預測的最小用氫成本在同一區塊上的位置，在對應實際生產條件范圍內尋找到最低用氫成本時的最優操作參數設置。