李元,章展鵬,徐鳴遠,陳丙珍,趙勁松
(化學工程聯合國家重點實驗室(清華大學),北京 100084)
報警系統建立的目的是防止化工生產過程中出現緊急停車、異常工況而引起經濟損失。但由于分布式控制系統(distributed control system,DCS)的出現,增設報警位點成本極低,這直接導致工廠報警數量急劇上升,無效報警過多,報警泛濫現象頻繁發生[1]。國際標準推薦操作員每10 min接受的平均報警數不超過1次,異常工況下10 min內的報警數不超過10次[2],而各領域報警系統的實際現狀與此標準差距很大。報警管理因此具有非常重要的地位。
隨著化工生產過程耦合度的不斷提高,相關報警也隨之增多。Rothenberg[3]指出,相關報警可分為序列報警與非序列報警,序列報警為某一位號發生報警所引起的固定報警序列,非序列報警為接連出現但延遲時間不固定的多個報警。對于相關報警的分析有利于報警系統更加簡潔、系統地展示報警信息,同時有助于操作員在此基礎上進行報警原因分析、異常工況預測等工作。
Izadi等[4]提出報警系統分析、設計、管理的標準與方法,指出合理的報警系統應在漏報率、誤報率、報警延遲等指標中做出平衡。Zhu等研究了重復報警的處置方法[5],并提出動態報警管理策略[6]。Choi等[7]提出了計算0、1序列之間的相似度,并利用凝聚層次聚類算法進行分組的聚類方法。接下來的研究工作均在此框架下展開。Chen等[8-10]提出報警相似度色彩圖(alarm similarity color map,ASCM)概念,認為它可以成為表征報警相似度的圖形工具。此外,有研究者對現有的相似度函數進行總結歸納,提出了相似度函數性能的評價標準[11-12],并指出可以使用互相關函數算法作為計算報警相似度的方法[13]。
現有的方法在處理序列報警時效果良好,但面對非規整報警時,傳統方法難以挖掘其相關性。這是因為現有方法利用整體平移0、1序列的方法引入延遲[11],卻沒有考慮到相似報警接連出現但延遲時間不固定的情況。其次,現有的方法均假定兩個報警位點的相關性對稱[12],即報警位點A對B的相關性等于B對A的相關性,但在實際生產中,會出現報警A必然導致報警B但報警B未必導致報警A的情況。這種非對稱的相關性是現有方法未考慮到的。
本文在傳統的基于凝聚層次聚類的報警相關性挖掘方法基礎上,提出了非規整報警相關性計算方法,彌補了傳統方法難以處理延遲時間變化與非對稱情況的不足。同時,為了能夠讓不同位點、不同裝置之間的相關性結果具有可比性,提出以概率形式表示報警位點之間的相關性大小的方法。最后,利用仿真案例與工廠真實生產案例驗證本方法的有效性。

圖1 典型的過程變量報警示意圖Fig.1 Typical sketch for process variable alarm
如圖1所示,在實際工業生產中,每一個位點設置有不同級別的報警限,當位點變量值超過報警限后,便會觸發相應級別的報警。為表征這種隨時間變化的報警序列,用“0”表示正常狀態,“1”表示超限(包括高報警限與低報警限)的報警狀態。以超出上限的情況為例,記過程變量序列為x(n),報警序列為xalarm(n),高報警限為xthreshold,則可以表示為

這樣的表征方法雖可以準確反映變量是否處于報警狀態,但在報警泛濫現象出現時,“1”的數量會急劇增多,導致需要挖掘的相關性被無效信息淹沒。因此,在實際應用中采用另一種方法,即僅將報警由正常變為超限的時刻記為“1”,其余均記為“0”

在得到各個位點的 0、1報警序列后,可利用相似度函數計算得到位點兩兩之間的相似度大小。接下來,利用凝聚層次聚類算法將位點之間的相似度聚類,得到若干組互相關聯的報警位點,并利用ASCM圖可視化地展示結果。
凝聚層次聚類算法是聚類算法的一種,由Johnson[14]首先提出,具體操作步驟如下:
(1)假設有N個對象,每個對象算作一個類C,得到類集合 {C1,C2,… ,CN},計算每對類兩兩之間的相似度,相似度計算方法可以自行定義,常見的是利用相似度函數進行計算(見1.3節);
(2)將相似度最大的兩個類Ci和Cj合并成一個新類CN+1;
(3)用新類CN+1替換掉Ci和Cj,得到新的類集合,且新集合中類的數量比舊集合中類的個數少1個;
(4)重新計算新集合中新類CN+1與其他類之間的相似度;
(5)重復步驟(2)~步驟(4),直至步驟(4)中計算得到的相似度低于設定閾值。最終得到的類集合即是聚類結果。
在步驟(4)中,需要定義新類CN+1與其他類之間的相似度如何計算。在本文中,計算兩個舊類Ci和Cj與其他類Ck之間的相似度,并將結果中的最大值作為新類CN+1與Ck的相似度。
將聚類后各對象之間的相似度記為維度為N的相似度矩陣,且同一類中的對象排列在對角線附近。以顏色的深淺表示相似度大小,則可以得到ASCM圖。以圖2為例,其中1~10表示10個對象,經過凝聚層次聚類后得到 4個類,并且C1= { 8},C2= { 1,6,7,4},C3= { 3,9},C4= { 10,2,5}。

圖2 ASCM示意圖Fig.2 Alarm similarity color map
Yang等[11]指出,針對報警系統性能最好的相似度函數為 Jaccard函數與 Sorgenfrei函數。其中,Jaccard函數的形式為

Sorgenfrei函數的形式為

式中,S為相似度大小,Na為序列a中出現的“1”的個數,Nb為序列b中出現的“1”的個數,C為兩個序列中共同出現的“1”的個數。在Na、Nb不變的情況下,C的值越大,S的值越大,也就是說兩個序列中同時出現的“1”越多,相似度越大。
為處理序列報警可能存在的整體延遲,可以對序列a進行時間長度為τ的整體平移,用平移得到的新序列aτ與b計算相似度S,并選取不同的τ多次計算以使S值最大,即

利用傳統的相似度函數可以計算有固定順序形式的序列報警之間的相關性。但在實際生產過程中,還會出現延遲時間不固定和非對稱的非規整報警,這類報警同樣具有相關性,卻是傳統的方法無法挖掘到的。
因此,本文給出如下步驟計算非規整報警序列的相關性:
(1)計算序列a中“1”的數量,記為N;
(2)設定窗口長度l;
(3)對于序列a中的每一個“1”,若在序列b對應位置前后各為l的長度內存在一個“1”,則記為一個“A1”,將序列a中“A1”的數量記為m;
(4)則令序列b對序列a的相關性 /Sm N= ;
(5)同樣的方法可計算序列 a對序列 b的相關性。
計算示例如圖3所示,令l= 1。對于序列a,位置P1處的“1”在前后l的范圍內不存在對應的“1”,而位置P2、P3處的“1”均可被記為“A1”,因此m=2 ,序列b對序列a的相關性S=2/3=0.66。

圖3 相關性計算示例Fig.3 Example for correlation calculation
由相似度函數得到的報警相關性大小受到報警數量的影響,此外,報警限的高低也會影響報警相關性大小,這導致不同過程、不同位號,利用相似度函數計算出來的相關性并不具備可比性,難以確定一個固定的標準來比較兩個報警序列的相關程度。
Yang等[11]指出,利用計算概率、不等式放縮與假設檢驗,可以得到傳統相關性函數的相關性閾值,進而確定兩個報警序列是否相關。但在使用非規整相關性計算方法后,序列a中一個位置的“1”是否記為“A1”不是單純取決于序列b中的對應位置,而是取決于 b中對應位置前后各為l的長度內是否存在“1”,因此無法套用 Yang的方法計算閾值。為此,本文提出如下方法衡量相關性的大小:
(1)對于待測的兩個序列a與b,記序列a中“A1”的個數為m;
(2)構造與待測序列中“1”的數量相同、但相互獨立的序列a′與b′,由于這樣的序列很多,因此記序列a′中存在的“A1”數量為隨機變量X;

其中,最核心的部分為計算概率分布P(X=x)。
令在序列a′中某一位置出現“1”的概率為ap′,序列b′中某一位置出現“1”的概率為bp′,并假設序列中各位置出現“1”的事件服從獨立同分布。定義事件I為:如果序列a′中特定位置O出現“1”,且序列b′中位置 O前后各長度為l的窗口中出現“1”,則認為出現事件I,同時序列a′中的“1”可記為一個“A1”。
則事件I出現的概率p為

求解概率分布P(X=x)可以等價于計算事件 I發生的次數。在實際生產過程中,pa′、pb′均遠小于1,即序列中“1”的個數極少,且l值往往較小,因此序列a′中多個位置的“1”對應于序列b′中同一個“1”的情況極為稀少,可以忽略。即認為事件I發生的概率服從獨立同分布,X服從二項分布

式中,n表示序列的長度。
為驗證式(7)的準確性,隨機生成 10000對相互獨立的0、1序列,序列長度均為2000,且有pa′=0.1,pb′= 0 .05,統計其中出現的 “A1”數量為k的序列對數量,并用此頻率值與式(7)計算得到的概率值比對,結果如圖4所示。

圖4 概率分布結果驗證Fig.4 Validation for result of probability distribution
可以看出,計算出的概率值與實際頻率基本吻合。
在實際計算式(7)的過程中,可以用“1”出現的頻率代替概率,即記序列長度為N,序列a′中“1”的個數為Na′,序列b′中“1”的個數為Nb′,則pa′=Na′/N,pb′=Nb′/N。
此外,可以仿照1.2節中介紹的ASCM圖,采用概率形式代替傳統計算方法得到的相關性,并將結果繪成色彩圖,這樣得到的圖像更加具有可比性。
假設i= 1 ,2,3,… ,1 00,而xA(i)是一個服從伯努利分布的隨機變量xp(i)∈ ( 0,1),且有P(xp(i)=0)= 0 .5,P(xp(i)= 1 )= 0 .5。對于t∈ [20i- 1 9,20i],構造一對長度為2000的相關序列x(t)和y(t)

式中,N(μ,σ) 表示均值為μ、標準差為σ的正態分布。
為序列x(t)設置報警限2.5[圖5 (a)],y(t)設置報警限3.0[圖5 (b)],利用式(2)的方法得到報警序列xa1(t)與ya1(t)[圖 5 (c)與圖 5 (d)],記為組 1#。將x(t)的報警限改為2.7,y(t)的報警限改為3.2,再次用式(2)的方法得到報警序列xa2(t)與ya2(t),記為組2#。此外,利用Matlab生成長度為2000的隨機0、1序列xa3(t)與ya3(t),顯然二者相互獨立,如圖6所示,記為組3#。分析這3組待測序列的相關性,得到結果如表1所示。

圖5 仿真案例相關序列Fig.5 Dependent sequence of simulation case

圖6 仿真案例獨立序列Fig.6 Independent sequence of simulation case

表1 仿真數據的相關性分析Table 1 Correlation analysis of simulation data
對比組1#與組2#的結果可以看出,在報警限改變后,非概率形式的相關性由0.29降為0.14。而利用概率形式表示的相關性分別為0.99與0.85,說明即使報警限改變,但計算結果依然可以說明待測序列相關。
對比組2#與組3#中非概率形式表示的相關性,前者為0.14,后者為0.15,這與組3#為獨立序列、組2#為相關序列的條件不符。這說明按照非概率形式,當報警限改變后,會出現相關序列計算出的相關性低于獨立序列的情況。而在概率形式表示下,組3#的相關性僅有0.46,相比于組1#與組2#明顯更低,是獨立序列。這樣的結果更符合實驗條件。
需要說明,在利用概率形式表示相關性大小時,得到的結果會整體偏大,如獨立序列組3#的相關性為0.46,遠大于非概率形式時的0.15。而事實上,組3#利用概率形式表示的相關性大小為0.46,意味著待測序列3ax中“A1”數量大于獨立序列3ax′中“A1”數量這一事件發生概率僅有46%,顯然不能說明待測序列3ax與3ay相關。
選取國內某化工廠聚合裝置作為分析對象,局部流程如圖7所示。以E釜為例,聚合釜R21E進料包括反應用水REC-WAT、引發劑INI,單體VCM經由泵P52A打入儲罐R52,再經由閥V28進入聚合釜。冷卻水COLD-WAT由泵P21E打入夾套,機封水SEAL-WAT進入釜底的攪拌機。聚合釜出料為氣相產物GAS與產品PRODUCT。

圖7 聚合工藝局部流程Fig.7 Scheme for part of polymerize process
選取該裝置中時間跨度為1個星期內報警最多的15個位號進行統計,利用傳統相關性計算方法得到ASCM圖,如圖8所示。做出使用非規整報警相關性計算方法并使用概率形式表示的色彩圖,如圖9所示。

圖8 傳統方法得到的ASCM圖Fig.8 Alarm similarity color map made by traditional method

圖9 利用非規整報警相關性計算方法且用概率形式表示的色彩圖Fig.9 Color map made by non-regular alarm correlation method in form of probability
可以看出,傳統的ASCM圖僅得到3組相關位點,分別為 FT-A21C與 C-A21C、FT-A21A與C-A21A、C-P31與FT-P31。這3組位點都為序列報警,且都是管路流量與該管路閥門閥位,其相關性顯而易見,因此指導意義有限。
而利用非規整報警相關性計算方法且用概率形式表示的色彩圖還可以得到以下3組相關位點。
(1)位點II-P52A對位點PT-R21E、PT-R21C、PT-R21D分別呈現出相關性。II-P52A為聚合釜單體進料泵的電流,而PT-R21E、PT-R21C、PT-R21D分別為3臺聚合釜的壓力。聚合釜單體進料增加會導致反應加劇,進而引發聚合釜的壓力升高,因此該相關性成立。而 PT-R21E、PT-R21C、PT-R21D對II-P52A的相關性并不成立,所以此組相關位點屬于非對稱的相關位點。
(2)位點PT-A21B與II-P21B互相呈現出相關性。II-P21B為聚合釜B的循環水泵電流,PT-A21B為聚合釜B的攪拌機封水壓力。聚合釜循環水泵加大功率后,會導致循環水流量上升,聚合釜溫度下降,聚合釜內以及與聚合釜聯通的機封水管路壓力均下降,因此該相關性成立。當釜內反應物的存量變化時,報警循環水泵功率過大與報警釜內壓力過低發生的間隔時間也會變化,所以此組相關位點屬于延遲時間變化的相關位點。
(3)不同釜之間攪拌機封水流量(如FT-A21A與FT-A21C)之間呈現出相關性。由于各分管路間的攪拌機封水流量由同一個總閥門控制,因此分管路間的攪拌機封水流量會同時波動,該相關性成立。當總管路與各分管路的流量變化時,某兩個分管路間攪拌機封水流量報警發生的間隔也會發生變化,所以此組相關位點屬于延遲時間變化的相關位點。
可以看出,利用非規整報警相關性計算方法且用概率形式表示的色彩圖可以得到更加全面、準確的結果。此外,上述15個位號中存在3組相關的非規整報警與3組相關的序列報警,相關的非規整報警組數占到全部相關組數的50%,這說明非規整報警相關性應該得到足夠的重視。
需要說明,聚合裝置的操作參數需要隨著產品牌號的不同而進行改變,而操作參數的改變必然會影響各位點的報警序列,進而影響到挖掘出的相關性結果。與此類似,當化工裝置的生產工況發生改變時,需要重新進行相關性挖掘。
為了進一步說明不同方法之間存在的差異,選取位點II-P52A與PT-R21D為待測序列組1,位點II-P21B與PT-A21B為組2,位點FT-A21C與C-21C為組3,C-A21B與II-P52A為組4,分別利用傳統的 Jaccard相關性函數、非規整報警相關性計算方法、概率形式表示的相關性大小進行分析。其中,組1為非對稱的相關位點,組2為延遲時間不固定的相關位點,組3為相關的序列報警位點,組4為相互獨立的對照組。在利用非規整報警相關性計算方法與概率形式時,僅分析每組兩個位點中前者對后者的相關性,如組2中即為II-P52A對PT-R21D的相關性。
此外,可以進一步計算組 1中 PT-R21D對II-P52A的相關性,得到利用概率形式表示的相關性大小為0.56。因此,PT-R21D對II-P52A的相關性不成立,這與組1為非對稱相關位點的條件一致。
由表2可以看出,非規整報警相關性計算方法與傳統相關性函數相比,可以更好地挖掘出延遲時間變化(組1)、非對稱(組2)的相關性位點。然而,僅利用非規整報警相關性計算方法而不使用概率形式表示時,組1與組4的結果十分接近,這與組1相關、組4獨立的條件不符。相比之下,利用概率形式表示的相關性結果更加合理。

表2 實際數據的相關性分析Table 2 Correlation analysis of industry data
本文在傳統的基于凝聚層次聚類的報警相關性發掘算法基礎上,對相關性計算方法進行改進,提出了非規整報警相關性計算方法,從而能夠更全面、更準確地挖掘化工過程中不同變量間的相關性。經由工業數據驗證,改進后的相關性發掘算法能夠挖掘出延遲時間變化、非對稱的相關位點,得到真實存在而傳統算法很難發掘到的相關性。
此外,本文提出了利用概率形式表示相關性大小的計算方法,使不同位號、不同裝置間的相關性大小具有可比性,這在3.2節的工業案例中得到了驗證。在報警相關性發掘算法中,一個重要的問題是當變量的報警限改變、報警數量大幅變化后,如何保證發掘相關性的算法結果依然成立。經過 3.1節中數據的驗證,利用概率形式表示相關性大小的方法可以有效解決此問題。
[1] Izadi I, Shah S L, Shook D S, Chen Tongwen. An introduction to alarm analysis and design//The 7th IFAC Symposium on Fault Detection, Supervision and Safety of Technical Processes [C].Barcelona, 2009.
[2] Engineering Equipment and Materials Users’ Association (EEMUA).EEMUA Publication 191:alarm systems—a guide to design,management and procurement [S]. London, 2007.
[3] Rothenberg D H. Alarm Management for Process Control:A Best-Practice Guide for Design, Implementation, and Use of Industrial Alarm System [M]. NewYork:Momentum Press, 2009.
[4] Izadi I, Shah S L, Shook D S, Kondaveeti S R, Chen Tongwen. A framework for optimal design of alarm systems//The 7th IFAC Symposium on Fault Detection, Supervision and Safety of Technical Processes [C]. Barcelona, 2009.
[5] Zhao Jinsong (趙勁松), Zhu Jianfeng (朱劍鋒). Data filtering based alarm processing strategy for repeating alarms [J].Journal of Tsinghua University:Science and Technology(清華大學學報:自然科學版), 2012, 52 (3):277-281.
[6] Zhu Jianfeng, Shu Yidan, Zhao Jinsong, Yang Fan. A dynamic alarm management strategy for chemical process transitions [J].Journal of Loss Prevention in the Process Industries, 2014, 30:207-218.
[7] Choi Seung-Seok, Cha Sung-Hyuk, Tappert C C. A survey of binary similarity and distance measures [J].Journal of Systemics,Cybernetics and Informatics, 2010, 8 (1):43-48.
[8] Kondaveetia S R, Izadi I, Shaha S L, Black T, Chen Tongwen.Graphical tools for routine assessment of industrial alarm systems [J].Computers & Chemical Engineering, 2012, 46:39-47.
[9] Ahmed K, Izadi I, Chen Tongwen, Joe D, Burton T. Similarity analysis of industrial alarm flood data [J].Automation Science and Engineering,IEEE Transactions on, 2013, 10 (2):452-457.
[10] Yang F, Shah S L, Xiao D, Chen T. Improved correlation analysis and visualization of industrial alarm data [J].ISA Transactions, 2012, 51(4):499-506.
[11] Yang Zijiang, Wang Jiandong, Chen Tongwen. Detection of correlated alarms based on similarity coefficients of binary data [J].Automation Science and Engineering, IEEE Transactions on, 2013, 10 (4):1014-1025.
[12] Faith D P. Asymmetric binary similarity measures [J].Oecologia,1983, 57 (3):287-290.
[13] Masaru N, Fumitaka H, Tsutomu T, Hirokazu N. Event correlation analysis for alarm system rationalization [J].Asia-Pacific Journal of Chemical Engineering, 2011, 6 (3):497-502.
[14] Johnson S C. Hierarchical clustering schemes [J].Psychometrika,1967, 32 (3):241-254.