韋咪娜, 袁德成, 莊亞文, 魏志斌

(沈陽化工大學 信息工程學院， 遼寧 沈陽 110142)

基于數(shù)據統(tǒng)計分析的因果關系建模研究

韋咪娜, 袁德成, 莊亞文, 魏志斌

(沈陽化工大學 信息工程學院， 遼寧 沈陽 110142)

針對復雜工業(yè)過程中變量多、變量間關系復雜、變量關系難以準確量化的特點，提出一種分析變量關聯(lián)關系及因果關系的數(shù)值方法，提高了因果關系分析的計算效率和精度.以TE化工過程模型為研究對象，針對TE流程中生產成本突增的問題，采集回路中控制偏差與生產成本的時間序列數(shù)據，建立其局部線性自回歸模型,分析多組控制偏差與成本間關聯(lián)性，運用數(shù)值算法計算了時域和頻域因果關系.結果表明:反應器壓力偏差是TE過程生產成本增加的主要原因，同時得到其他變量間的關聯(lián)性.通過仿真實例驗證了方法的有效性.

因果關系模型； 復雜工業(yè)過程； 時間序列數(shù)據

基于數(shù)據驅動的因果關系建模在大型復雜工業(yè)過程的分析和設計中發(fā)揮著越來越重要的作用，尤其是在解決系統(tǒng)異常根本原因分析、異常報警管理、全流程控制系統(tǒng)設計方面具有重要應用價值[1].在全流程復雜系統(tǒng)，如現(xiàn)代工業(yè)過程、生物系統(tǒng)、社交網絡，各變量間互相聯(lián)系，系統(tǒng)行為取決于每對變量間的相互關系及每一變量的局部動態(tài)特性，因此，識別這些關聯(lián)關系，即關聯(lián)性與因果關系對于分析系統(tǒng)的影響機制、結構特征和整體動態(tài)特性非常重要.連通性和因果關系獲取方法有知識驅動和數(shù)據驅動兩種形式：知識驅動作為一種依據過程知識定性建模的方法，由于對過程知識的完整性要求十分高，因此，并非總是可用；而現(xiàn)代化工業(yè)企業(yè)每天運行在DCS數(shù)據庫中，積累了大量數(shù)據，為建模及分析提供了寶貴資源，因此，建立基于數(shù)據驅動的拓撲結構，構建因果關系網絡更加實用可行，這已成為學術界和工業(yè)界關注的重點.

針對過程變量，探測因果關系的數(shù)據驅動方法分為三類：滯后模型法(包括格蘭杰因果關系、傳遞熵)、條件獨立法(貝葉斯學習)、高階統(tǒng)計量法.Duan P等運用傳遞熵測定有向因果關系[2]，雖然傳遞熵法無模型、無線性限制，但存在過渡依賴PDFs(工藝物料平衡圖)，增大了計算負擔，時間滯后無法估計，且易受非平穩(wěn)噪聲影響.Baucer M等運用互相關分析識別廠級干擾傳播路徑，通過建立CCF(互相關函數(shù))彌補了傳遞熵法無法估計時間滯后的不足[3]，其計算方便，但忽視了時間序列趨勢，所得時間滯后僅為估計值，計算精度仍有待改進.Cowell R G等通過貝葉斯學習建立隨機網絡及專家系統(tǒng)[4]，盡管循環(huán)因果關系算法的開發(fā)彌補了傳統(tǒng)貝葉斯網絡無環(huán)圖在獲取系統(tǒng)動態(tài)特征中的不足，但模型的建立需要所有模式下的數(shù)據，甚至包括所有故障模式下的數(shù)據，成為建立貝葉斯網絡的最大限制.

本文通過建立時間序列數(shù)據的多變量的向量自回歸模型(VAR)進行格蘭杰因果分析(Granger Causality Analysis)，在時域分析的基礎上，對因果影響進行譜分解，對多元自回歸過程的分解使得因果關系的譜表示更加明確有效，采用數(shù)值算法建立VAR模型，精確估計了模型的階與滯后系數(shù)，在時間序列數(shù)據長度足夠的情況下可檢測到變量間微弱相關關系，并在可接受計算負擔內提高了計算效率和精度.實驗仿真證明了本文運用數(shù)值途徑實現(xiàn)格蘭杰因果分析方法的有效性.

1 數(shù)據獲取與處理

美國Tennesses Eastman化工公司在實踐基礎上研究出了一個可進行控制方案設計、控制策略研究的仿實際過程的工業(yè)模擬模型，即TE化工過程模型.TE過程模型變量多，變量關系復雜，且連續(xù)系統(tǒng)中各單元間有較強的相互作用，是一個典型的復雜多約束非線性過程[5].本文選取TE過程9個控制回路的控制偏差與生產成本共10個變量為研究對象，在TE過程的simulink模型中加入隨機干擾，進行3次試驗，運行TE仿真模型采樣，得到維數(shù)為3×1 000×10的數(shù)據，即10個過程變量的時間序列.

基于數(shù)據分析的因果關系建模對序列的穩(wěn)定性非常敏感，非平穩(wěn)時間序列的因果關系檢驗會得到偽回歸，導致產生虛假因果關系[6].因此，在進行因果關系建模前，需對采集到的數(shù)據進行平穩(wěn)性檢驗.對于時間序列Xt可通過下式進行平穩(wěn)性檢驗：

(1)

ΔXt為t時刻序列數(shù)據變化量，δ為時間序列系數(shù)，εt表示噪聲序列.原假設和備擇假設為H0:δ≥0;H1:δ<0;通過t統(tǒng)計檢驗量進行檢驗，當t統(tǒng)計量小于τ，拒絕原假設H0，Xt序列平穩(wěn)；若t統(tǒng)計量大于τ，則接受原假設，認為Xt存在單位根，為非平穩(wěn)序列.如果序列為非平穩(wěn)序列，接下來繼續(xù)檢測ΔXt的平穩(wěn)性，直至檢驗結論表明差分后的序列為平穩(wěn)序列為止.

2 因果模型建立

2.1 VAR模型的建立

建立向量自回歸模型(VAR)是目前研究變量間格蘭杰因果關系的主要方法.針對兩個隨機過程Xt和Yt，可分別建立它們的VAR模型：

(2)

(3)

單變量向量自回歸模型反映了隨機過程當前值與過去值之間的關系，方差Σ1或Γ1表明這種預報關系的正確程度.也可建立兩個隨機過程間的聯(lián)合回歸模型，反映兩隨機過程間的關聯(lián)性.

(4)

(5)

其中噪聲誤差項在時間上不相關，在同一時刻的自相關矩陣定義為:

(6)

式中:Σ2=var(ε2t),Γ2=var(η2t),如果Xt與Yt不相關，則{b2j}和{c2j}等于零，γ2=cov(ε2t,η2t)等于零，Γ2=Γ1.

在復雜系統(tǒng)中，并非總是僅存在兩變量關系，無約束回歸將導致虛假因果關系的產生.當隨機過程Yt與Xt間不存在直接因果關系，但Xt、Yt都與隨機過程Zt存在滯后相關關系時，無約束回歸將會錯誤地預報Xt、Yt間的因果關系[7].在此情況下，建立3變量間的有約束回歸模型將會有效消除因果推斷中一些變量可能造成的混雜影響.對于3個隨機過程Xt、Yt和Zt，先建立Xt和Zt之間的VAR模型：

(7)

(8)

式中協(xié)方差矩陣為：

(9)

進一步加入第3個隨機過程Yt，建立3個過程間的VAR模型：

ε4t, var(ε4t)=Σxx

(10)

η4t, var(η4t)=Σyy

(11)

χ4t, var(χ4t)=Σzz

(12)

其中，協(xié)方差矩陣

(13)

當Yt與Xt間的因果關系可以完全由Zt間接表達時，{b4j}為零矩陣，Σ3=Σxx，則在Zt約束下，Yt對Xt的因果關系為零，加入隨機過程Yt的過去值對Xt的預測效果沒有改善.另一方面，如果Yt對Xt存在直接影響，加入Yt的過去值將會更好地預測Xt值，Σ3>Σxx，存在在Zt約束下Yt對Xt的因果關系[8].

對于所建立的VAR模型，運用杜賓-瓦特森(Durbin-Watson)統(tǒng)計量檢測回歸分析中的殘差項是否存在自相關.根據模型估計數(shù)據預測時間序列趨勢，與實際數(shù)據對比，檢驗模型一致性.

2.2 時域格蘭杰因果關系分析

兩個隨機過程Xt和Yt之間的總體相關性為：

(14)

Σ表示封閉(有限項)矩陣的行列式，Σ1或Γ1分別為兩隨機過程方差.當兩個時間序列Xt和Yt不相關時，它們的總體互相關性FX,Y=0成立；反之，則有FX,Y>0.參照式(2)、(4),如果Σ2小于Σ1，則可以判斷加入Yt的過去項后Xt預報更加精確，也可以認為Yt對Xt因果影響.它們之間的這類因果影響關系為：

(15)

如果FX→Y等于零，從Yt到Xt沒有因果關聯(lián)；如果FX→Y大于零，則存在相互影響.類似地，Xt到Yt的因果關聯(lián)如下：

(16)

Xt到Yt的間接相關性：

(17)

如果γ2=cov(ε2t,η2t)等于零，則FX,Y等于零；反之，F(xiàn)X,Y大于零.根據以上定義，兩變量之間的因果關系可以統(tǒng)一在一個公式中：

FX,Y=FX→Y+FY→X+FX·Y

(18)

引入第3個變量后，可能存在直接和間接的關聯(lián)，如圖1所示.

圖1 三變量因果關系

經過Zt傳遞的、從Yt到Xt的因果關聯(lián)為：

(19)

2.3 頻域格蘭杰因果關系分析

格蘭杰因果關系具有將因果關系進行相應頻域分解的突出優(yōu)勢，能夠清晰反映兩個隨機過程產生關聯(lián)的頻率[10].頻域分解得到的譜因果分析整合了時域的因果關系，可認為是所有頻率譜因果關系的平均值.隨機過程Xt、Yt頻域分解的譜因果關系為:

fY→X(λ)≡

(20)

Sxx(λ)=Hxx(λ)Σxx(λ)Hxx(λ)*+

Hxy(λ)Σyy(λ)Hxy(λ)*

(21)

偏協(xié)方差矩陣：

Σy|x≡Σyy-ΣyxΣxx-1Σxy

進一步針對3個隨機過程Xt、Yt、Zt建立頻域內受約束回歸模型.先建立Xt、Zt間VAR模型：

(22)

(23)

X*、Z*為白噪聲殘差，其與隨機過程Xt的轉換如下:

(24)

則頻域內3個隨機過程的有約束的譜因果關系為

fY→XZ(λ)≡fY⊕Z*→X*(λ)

(25)

3 實驗仿真

3.1 VAR模型參數(shù)

VAR模型參數(shù)：模型階數(shù)為7階，滯后系數(shù)5 975，譜半徑 0.996 922，相對誤差3.279 9×10-11，時域與頻域因果關系最大相對誤差9.68×10-12<1×10-5.模型一致性100 %.杜賓-瓦特森統(tǒng)計檢測(D-W檢驗)結果見表1.

表1 VAR模型D-W統(tǒng)計檢驗

可以看出變量殘差在0.05的顯著性水平上不存在自相關.修正后的R2基本大于90 %，模型擬合良好.

3.2 因果關系仿真結果

格蘭杰因果關系及顯著性分析如圖2、圖3所示.

圖2 時域格蘭杰因果關系

圖3 格蘭杰因果關系顯著性檢驗

從圖2可以看出：在時域內建立10維時間序列數(shù)據的多約束自回歸模型，變量2→10存在明顯格蘭杰因果關系，其次變量7→10、7→6、7→2 等也存在單向關聯(lián)性.相應地，這些因果關系或關聯(lián)關系在網格圖中顯示P值檢驗極小，顯著性檢驗結果顯著，驗證了格蘭杰因果關系分析結果的可靠性.在圖3中，進一步對回歸模型中每對關聯(lián)關系顯著性進行了理論檢驗與置換檢驗.置換檢驗作為基于計算機的模擬方法，不依賴原始數(shù)據分布，以兩樣本均數(shù)之差作為統(tǒng)計量，充分利用樣本數(shù)據信息，提高了檢驗效能并再次驗證了所得因果關系.

在時域因果關系分析的基礎上，編寫格蘭杰因果分析算法程序,將隨機過程進行譜分解，變換時域中VAR模型參數(shù)為頻域中新變量，并根據式(25)計算譜因果關系.分析10組隨機過程變量在GC(Granger Causality)因果關系最大的頻率點上的格蘭杰因果關系連接，發(fā)現(xiàn)10組時間序列數(shù)據兩兩配對因果關系呈現(xiàn)出明顯峰值，在10 Hz頻率點附近變量2→10、7→10、7→6之間的單向因果關系顯著，仿真結果如圖4所示.從圖4與圖5可以看出：變量2(反應器壓力偏差)與7(汽提塔出流量控制偏差)對10有直接作用，變量3與2、7與2等存在微弱的相關關系.

圖4 譜因果關系

圖5 零分布的K-S檢驗圖

變量間關聯(lián)關系見圖6.置換檢驗結果表明所得譜因果關系是顯著的，進一步證實了反應器壓力偏差、汽提塔出流量控制偏差這兩個變量與生產成本間的單向因果關系.TE過程的經濟消耗主要由原料消耗、放空氣體中殘留產品和原料及處理F產品的消耗組成，反應器壓力控制偏差可直接影響物料反應速率與轉化率，未完全轉化原料的流失會增加經濟消耗，因此，反應器壓力控制偏差與生產耗費間的因果關系不違背實際生產情況，這也與以上時域因果分析的結果一致.

圖6 變量間關聯(lián)關系

3.3 仿真驗證

運行仿真系統(tǒng)，驗證反應器壓力的測量值與生產成本之間因果關系，如圖7、圖8所示.

圖7 反應器壓力圖

圖8 TE過程生產成本圖

其中圖7中橫坐標為樣本時間序列數(shù)據，采樣單位為s，縱坐標表示反應器壓力值，單位為kPa；圖8中橫坐標為樣本時間序列數(shù)據，采樣單位為s，縱坐標表示生產成本值，單位為萬元/最大生產率.

圖7與圖8中曲線1和曲線2分別為第一次、第二次實驗測量值，兩次實驗中反應器壓力設定值均為2 800 kPa.當反應器壓力偏離設定值，生產成本隨之增加，并且偏差越大，成本增加越多.反應器壓力偏差在接近第100個采樣單位處達到最大，生產成本在第100個采樣單位后亦達到最大值.觀察圖7、圖8可得反應器壓力偏差與生產成本的強因果關系，并存在約0.5個單位時間的滯后，驗證了該方法的有效性.

4 結 論

建立了時間序列數(shù)據的VAR模型，運用數(shù)值算法實現(xiàn)對多變量關聯(lián)性及因果關系的分析.得出時域和頻域中因果關系，兩者最大絕對誤差9.68×10-12，分析結果基本一致.運用數(shù)值解法計算變量因果關系，提高了計算效率和計算精度，可以解決統(tǒng)計不準確和計算效率低的常見問題.仿真驗證表明：存在反應器壓力偏差與生產成本間的因果關系，證明了該方法的正確性.

[1] YANG F,SHAH S,XIAO D Y.Signed Directed Graph based Modeling and its Validation from Process Knowledge and Process Data[J].Int J App Math Comput Sci,2012,22(1):41-53.

[2] DUAN P,YANG F,CHEN T，et al.Direct Causality Detection via the Transfer Entropy Approach[J].IEEE Trans Control Syst Technology,2013,21(6):2052-2066.

[3] BAUER M,THORNHILL N F.A Practical Method for Identifying the Propagation Path of Plant-wide Disturbances[J].J Process Control,2008,18(7):707-719.

[4] COWELL R G,DAWID A P,LAURITZEN S L.Probabilistic Networks and Expert Systems[J].Springer-verlag,1999,39(20):2548-2558.

[5] 袁德成.過程控制工程[M].北京：機械工業(yè)出版社，2013：13-21.

[6] OYELEYE O O,KRAMER M A.Qualitative Simulation of Chemical Process Systems:Steady-State Analysis[J].AICHE J,1988,34(9):1441-1454.

[7] YANG F,XIAO D Y.Review of SDG Modeling and its Application[J].Control Theory & Applications,2005,22(5):767-774.

[8] SETH A K.A Matlab Toolbox for Granger Causal Connectivity Analysis[J].J Neurosci Methods,2010,186(2):262-273.

[9] YIM S Y,ANANTHAKUMAR H G,BENABBAS L,et al.Using Process Topology in Plant-wide Control Loop Performance Assessment[J].Comput Chem Eng,2006,31(2):86-99.

[10] GRANGER C W J.Investigating Causal Relations by Econometric Models and Cross-spectral Methods[J].Econometric,1969,37(3):424-438.

Abstract: A great number of process variables existed in the complex industrial process.The relationship between these variables were complicated and very difficult to quantify accurately.A numerical method of analyzing the variable relationship and causality was put forward aimed at improving the calculation efficiency and precision.Focus on TE process model,a great deal of time series data of control deviation and production cost was collected to establish their local linear auto regression model,in order to interpret the soaring production cost.Connectivity among these variables was discussed and some numerical algorithm was adopted to compute causality relationship in both time domain and frequency domain.The result showed that the reactor pressure deviation is the main cause of the raise of the cost of TE production process.Simultaneously,the connectivity and correlation of other variables were captured.The simulation validated this method was correct and effective.

Keywords: causal model; complex industrial process; time series data

StudyontheCausalModelingBasedonStatisticalAnalysisofData

WEI Mi-na， YUAN De-cheng， ZHUANG Ya-wen， WEI Zhi-bin

(Shenyang University of Chemical Technology, Shenyang 110142， China)

10.3969/j.issn.2095-2198.2017.03.015

TP1

A

2015-03-24

韋咪娜(1989-)，女，河南洛陽人，碩士研究生在讀，主要從事數(shù)據驅動、統(tǒng)計分析等方面的研究.

袁德成(1960-)，男，內蒙古阿拉善人，教授，博士，主要從事化工過程建模、控制與優(yōu)化的研究.

2095-2198(2017)03-0266-07