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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

2 因果模型建立

2.1 VAR模型的建立

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

(2)

(3)

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

(4)

(5)

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

(6)

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

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

(7)

(8)

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

(9)

進(jìn)一步加入第3個隨機(jī)過程Yt，建立3個過程間的VAR模型：

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

(10)

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

(11)

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

(12)

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

(13)

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

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

2.2 時域格蘭杰因果關(guān)系分析

兩個隨機(jī)過程Xt和Yt之間的總體相關(guān)性為：

(14)

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

(15)

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

(16)

Xt到Y(jié)t的間接相關(guān)性：

(17)

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

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

(18)

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

圖1 三變量因果關(guān)系

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

(19)

2.3 頻域格蘭杰因果關(guān)系分析

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

fY→X(λ)≡

(20)

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

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

(21)

偏協(xié)方差矩陣：

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

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

(22)

(23)

X*、Z*為白噪聲殘差，其與隨機(jī)過程Xt的轉(zhuǎn)換如下:

(24)

則頻域內(nèi)3個隨機(jī)過程的有約束的譜因果關(guān)系為

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，時域與頻域因果關(guān)系最大相對誤差9.68×10-12<1×10-5.模型一致性100 %.杜賓-瓦特森統(tǒng)計檢測(D-W檢驗)結(jié)果見表1.

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

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

3.2 因果關(guān)系仿真結(jié)果

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

圖2 時域格蘭杰因果關(guān)系

圖3 格蘭杰因果關(guān)系顯著性檢驗

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

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

圖4 譜因果關(guān)系

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

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

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

3.3 仿真驗證

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

圖7 反應(yīng)器壓力圖

圖8 TE過程生產(chǎn)成本圖

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

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

4 結(jié) 論

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

[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].北京：機(jī)械工業(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ù)據(jù)驅(qū)動、統(tǒng)計分析等方面的研究.

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

2095-2198(2017)03-0266-07