基于慢特征分析的分布式動態工業過程運行狀態評價

鐘林生 常玉清 王福利,2 高世紅

隨著科技水平的快速發展,企業對生產過程的精細化管理提出了更高要求.如何在保證生產安全和產品質量的前提下進一步提高綜合經濟指標,成為企業在當今激烈競爭的市場環境中存活的關鍵.為實現上述目標,企業不僅要求生產過程處于正常的狀態,還渴望追求過程運行在優狀態[1].在企業投產初期,工程師會利用過程知識和優化技術,使得過程處于優運行狀態.但是,由于外部擾動、過程參數的漂移和現場人員的不當操作等原因,過程會逐漸偏離最初設定的優運行區域[2].如果不能及時地檢測出過程運行性能的退化,輕則降低企業的經濟效益,重則演化為故障,威脅企業生產人員的安全.為了應對以上問題,過程運行狀態評價和過程監測技術應運而生.運行狀態評價是評估工業過程運行狀態優劣等級的方法,其隸屬于廣義上的過程監測范疇,但在實現目標方面又有顯著不同[3].具體而言,與傳統過程監測技術僅將過程劃分為正常與故障兩種狀態相比,運行狀態評價方法進一步將正常運行的生產過程細分為多個性能等級(例如優和非優).運行狀態評價可以看作是能夠滿足更高企業要求的擴展型過程監測方法.因此,研究及時、準確和全面的工業過程運行狀態評價方法,對保障企業生產安全和提高企業綜合經濟指標,具有重要的理論價值和實際意義.

在過去的20 年中,運行狀態評價獲得了學術界和工業界的廣泛關注.為了處理工業過程的多模態特性,Ye 等[4]采用高斯混合模型識別過程所處的模態,并開發了一種將優性評價指標劃分為不同等級的分類方法.考慮到并非所有過程變量都會顯著地影響過程的運行狀態,Fan 等[5]首先采用皮爾遜相關系數選出與評價指標高度相關的變量,然后利用最小二乘支持向量機模型實現過程運行狀態的劃分.針對新過程訓練數據不足問題,Zou 等[6]采用交叉域特征遷移算法,將舊過程的通用信息遷移到新過程中,從而提高新過程的評價模型的精度.為了同時評價過程穩定模態和過渡模態的運行狀態,Liu 等[7]基于綜合經濟指標制定評價策略,對穩定模態通過高斯過程回歸預測綜合經濟指標,對過渡模態采用歷史數據庫中相似過渡模態的綜合經濟指標的加權平均值作為當前過渡模態的評價指標.為了更好地處理過程數據中存在的不確定性信息,Wang 等[8]結合特征選擇、即時學習和概率主成分回歸,建立評價指標的預測模型,并通過概率模糊推理將多個評價指標轉換為最終的評價結果.雖然文獻[1-8]對過程運行狀態評價進行了廣泛探索并取得了巨大成就,但是上述研究成果普遍基于靜態模型,忽略了過程的動態特性.當過程數據中包含動態信息時,直接將靜態模型應用到動態數據是一種線性靜態近似,并不能準確地揭示變量之間的關系[9].在廣泛應用閉環控制的工業過程中,變量的時序相關性是普遍存在的[3].為挖掘過程數據的動態特性,學者們相繼提出了典型變量分析、動態獨立成分分析、動態主成分分析(Dynamic principal analysis,DPCA) 和動態慢特征分析(Dynamic slow feature analysis,DSFA)等方法[9-12].Sedghi等[13]采用動態主成分回歸技術預測過渡模態的綜合經濟指標,改善了文獻[7]的運行狀態評價性能.Zou 等[14]將典型變量分析和慢特征分析(Slow feature analysis,SFA)進行分層組合,提取過程更深層次的特征,從而更加細致地刻畫運行狀態的等級.雖然上述方法[9-14]對一般的動態過程很有效,但是它們都是集中式模型,在應用于規模較大的工業過程時,很可能會因為忽略了過程的局部信息而降低模型評價非優狀態的靈敏性.需要注意的是,上述基于分類[4-6]或回歸技術[7-8,13]運行狀態評價方法的應用前提是能獲得足夠多優狀態和非優狀態的訓練數據.但在實際生產過程中,任意一個或多個變量發生不同程度的偏移都有可能導致過程處于非優狀態.由于經濟和技術條件限制,在企業投產初期,幾乎不可能獲得所有非優狀態的完整數據集.針對非優狀態模式多樣但可獲得的非優狀態數據量相對較少的情況,用于表示優運行區域的閉合決策邊界往往會比表示優與非優狀態的分隔決策邊界更具可靠性.以多元統計分析技術[15]和一類分類技術[16]為代表的數據描述方法為上述情況提供了很好的解決方案.

針對現代工業過程規模大、流程長和工序多的特點[17],分布式模型得到了越來越多學者們的關注.分布式建模策略是經過過程分解、子模型建立和子模型結果融合而得到,因而更適用于規模較大的工業過程[18].Liu 等[19]提出一種基于全潛結構投影模型的分布式運行狀態評價方法.為了應對工業過程大數據問題,Zhu 等[20]基于MapReduce 框架,設計一種分布式并行主成分分析方法.考慮到字典學習能夠有效減少高維數據的計算和存儲負擔,Huang等[21]提出一種基于字典學習的分布式過程監測方法.然而,上述方法[18-21]均根據過程拓撲或機理知識實現過程的分解,有可能無法將過程變量劃分為概念上有意義的重疊或不相交的子塊[22].此外,當面對大規模和高復雜度的工業過程時,很難保證過程知識的完備性.為了解決這個問題,近年來涌現出一些數據驅動的過程分解方法[23],即根據過程變量的統計特性將其劃分為多個子塊.Ge 等[24]根據主成分的不同投影方向構造子塊,將原始特征空間自動地劃分為若干個子特征空間.為了處理過程的非線性,Jiang 等[25]利用互信息將相似性高的變量劃分為同一子塊,并對每一子塊建立核主成分分析模型,從而有效地挖掘出過程數據中的非線性信息和局部信息.顯然,以上過程分解方法[18-21,24-25]并未充分利用過程的動態信息.為了解決該問題,Tong等[22]基于皮爾遜相關系數,提出一種動態特征選擇方法,但其子塊個數等于過程變量數,這在一定程度上過度強調過程的局部信息,很可能會導致模型的誤報率升高.Zhong 等[26]在文獻[22]基礎上,采用最小冗余最大相關性(Minimal redundancy maximal relevance,mRMR)實現過程的分解,但mRMR 本質上是一種非線性特征選擇方法,并不能充分利用過程數據的動態信息.

以上分析表明,充分利用動態信息的過程分解方法還需要進一步研究.此外,過程分解方法需要與相應的子模型的建模原理密切相關,才有可能顯著地提高分布式評價模型的性能.以此為動機,本文提出一種分布式動態過程運行狀態評價方法.首先,結合動態時間規整(Dynamic time warping,DTW)[27]和K-medoids 聚類算法[28],實現過程子塊的劃分;然后,針對每一變量子塊建立相應的DSFA模型;最后,利用貝葉斯推理獲得全局的綜合評價指標,從而實現大規模工業過程的運行狀態評價.為方便后續闡述,將本文方法命名為分布式動態慢特征分析(Distributed dynamic slow feature analysis,DDSFA).本文工作的主要貢獻如下: 1)結合動態時間規整和K-medoids 聚類算法對過程進行分解,將波形相近的過程變量劃分到相同的子塊,提高了DSFA 算法對過程動/靜態信息的表征能力;2)與傳統動態模型為所有變量選擇單一時滯系數不同的是,本文設計一種新穎的時滯系數確定方法,使得不同子塊可以具有不同的時滯系數,增加了子模型的多樣性.

1 預備知識

1.1 動態時間規整

動態時間規整是一種常用的度量2 個時間序列相似度方法,主要是形狀相似度[27].通過對2 個時間序列進行拉伸或壓縮,盡可能地消除序列在時間軸上的扭曲,從而找到序列間相似度最大的匹配路徑.給定如下2 個時間序列:

式中,I和J為2 個時間序列的長度.需要注意的是,2 個時間序列的長度可以是不相等的,但2 個序列中元素的維度必須相同.為了便于理解,以下僅考慮一維時間序列.DTW 的目標是尋找一個規整路徑F=[f1,f2,···,fK],使得序列a和b的累積距離最小,其中fk=(ik,jk) 表示序列a的第ik個元素與序列b的第jk個元素相匹配.通常采用歐氏距離的平方來度量2 個序列元素間的差異,即:

DTW 的優化目標函數可表示為:

為了盡可能保留序列在時間軸上的本質結構并確保能夠找到最優匹配路徑,在尋找規整路徑過程中,還需滿足以下約束條件:

1)單調連續性條件.單調連續性條件是為了保證所獲得的規整路徑沒有跳過任何元素,也沒有出現交叉情況:

2)邊界條件.在規整路徑F中,2 個序列的開始和結束元素必須相互匹配:

3)窗口調整條件.增加調整窗口長度條件,一是考慮到2 個時間序列在時間軸上的扭曲通常不會太大,二是為了縮短路徑搜索的時間:

由式(4)～ (7)定義的優化問題可以采用動態規劃算法獲得最優解.d(a,b) 取值越小,代表序列a和b的相似度越大;反之亦然.

1.2 慢特征分析

慢特征分析是一種常用的降維技術,其前提假設是變化緩慢的特征比變化快速的特征更能代表數據核心信息[12].該觀點已被證實與生物視覺特性一致[29].SFA 的優化目標是尋找一個映射函數使得轉換后的特征的變化盡可能慢.假設存在一個m維的輸入信號x(t)=[x1(t),x2(t),···,xm(t)]T,在線性情況下,慢特征s(t) 可表示為:

SFA 的優化目標函數定義如下:

式中,〈·〉t表示時間平均運算符,為慢特征sj對時間的一階導數.需要注意的是,每個慢特征sj都被約束為零均值和單位方差,從而避免常數解.此外,不同慢特征之間的去相關,使得不同慢特征所攜帶的信息不同.

由式(9)的約束條件,可得:

由式(11)、式(12),可得:

與DPCA 類似,可以直接將上述SFA 算法應用于原始輸入變量的時滯拓展矩陣,來獲得其動態版本DSFA.時滯拓展矩陣進一步考慮了每個變量的l個滯后采樣時刻,并按如下方式進行疊加[11]:

2 分布式運行狀態評價方案

通過對DTW 和DSFA 方法的分析表明,DSFA算法的目標是從過程數據中提取出具有較小變化速度的潛變量,即潛變量波動盡可能平緩.顯然,潛變量變化速度與其自身波形相關.反向思考,波形相近的過程變量很可能是由相同潛變量所驅動的.因而,如果將波形相近的過程變量劃分到相同變量子塊,則有望提高所提取潛變量的準確性.考慮到,DTW 算法盡可能消除了2 個序列間在時間軸上的扭曲,是度量2 個時間序列波形相似性的常用方法.因此,本節創新性地將DTW 和DSFA 進行組合,設計一種分布式運行狀態評價方案,主要步驟包括基于DTW 和K-medoids 聚類算法的過程分解、基于DSFA 的運行評價子模型的建立和基于貝葉斯推理的子塊評價結果融合.

2.1 基于DTW 和K-medoids 聚類算法的過程分解

傳統的靜態過程分解算法通常假定訓練集中樣本是獨立同分布的.然而,工業過程的動態特性會導致過程變量存在時序相關性(即自相關和互相關)[30].現有過程分解算法較少考慮過程的動態特性,因此還需要進一步研究利用過程數據動態信息的過程分解算法.給定一個標準化的過程運行在優狀態的訓練集X=[x1,x2,···,xm]T∈Rm×n,其中xi∈R1×n.m和n分別代表變量數和樣本數,則過程變量間的DTW 距離矩陣D ∈Rm×m可構造為:

式中,Dij=d(xi,xj),表示變量xi和xj的DTW距離.考慮到,K-medoids 選取的質心是真實存在的過程變量,能有效避免虛擬質心對簇內過程變量間波形相似度的影響.基于距離矩陣D,本文采用K-medoids 聚類算法,將過程變量劃分為G(G

1)從過程變量x1,x2,···,xm中隨機選出G個變量作為簇的中心O1,O2,···,OG.

2)重復下面過程直到算法收斂.

a) 確定變量的所屬簇C={C1,C2,···,CG}.通過式(17),將每個變量分配到與簇中心的DTW距離最小的簇:

b)更新簇的中心.確定所有變量所屬的簇后,通過式(18)重新選擇每個簇的中心.其中 card(·) 是計算集合元素個數的運算符:

2.2 基于DSFA 的運行狀態評價子模型的建立

第2.1 節已將過程變量劃分為G個子塊X=[X1,X2,···,XG]T,針對每個子塊Xg ∈Rmg×n,建立其對應的DSFA 模型.動態建模方法的關鍵步驟是時滯系數l的確定.大部分現有的動態建模方法都為所有變量選擇了相同的時滯系數.然而,這種做法一方面,人為增大了建模數據的維度[31],從而導致模型的計算復雜度增大;另一方面,規模較大的工業過程通常包含多個子工序,每個子工序都具有不同的運行機制和控制策略,因此,過程變量一般不太可能都具有相同的時滯結構[32].綜合以上考慮,本文提出一種針對每個子塊變量特性時滯系數lg的確定方法.由第2.1 節可知,DTW 的目標是尋找使得序列間相似度最大的匹配路徑.因而,在獲得變量間的DTW 距離時,還同時獲得了變量間的最佳匹配路徑.具體地,假設同一子塊內變量xi和xj的 最佳匹配路徑為Fij=[f1,f2,···,fK],其中fk=(ik,jk) 表示變量xi ∈Xg的第ik個元素與變量xj ∈Xg的第jk個元素相匹配.在實際工業過程中,過程數據通常是由等間隔采樣獲得.此外,訓練集中不同變量的樣本數相同.綜合以上分析,變量xi和xj是等長的時間序列.因而,其時滯系數可以粗略估計為 |ik-jk|.但由于過程噪聲的影響,而且在實際過程中,高度相似的變量相對較少.因而,僅依靠最佳匹配路徑Fij中的一個元素確定2 個變量的時滯系數是不可靠的.鑒于此,本文提出一種更具一般性的時滯系數確定方法:

式中,m ode{·} 表示取眾數運算符,|?Fij|={|ikjk||(ik,jk)∈Fij},δ是一個非負整數.增加δ項,一方面,是為了更全面地提取過程數據蘊含的動/靜態信息;另一方面,是為了便于在實際應用中對式(19)進行細調.

通過式(19)確定每一子塊的時滯系數lg后,依據式(15)獲得對應子塊Xg ∈Rmg×n的時滯拓展矩陣Xg(lg)∈Rθ×N,其中,θ=mg(lg+1).利用SFA算法提取Xg(lg) 的慢特征,并根據特征的變化速度將慢特征劃分為以下2 個部分:

式中,sg,d ∈Rr×N是變化較慢的特征,包含了系統的核心信息;sg,e ∈R(θ-r)×N是變化較快的特征.通過式(22)的分位數方法確定r值[12]:

式中,Qq{〈〉t} 表示集合 {〈〉t} 的q分位數.簡言之,主導子空間sg,d移除了比所有原變量變化速度q分位數大的特征.

為指示主導子空間和剩余子空間中過程穩態分布p(x) 的變化,構造以下2 個統計量:

相應地,系統對過程暫態分布p() 的變化可通過以下2 個統計量指示:

2.3 基于貝葉斯推理的子塊評價結果融合

給定一個標準化的在線樣本xnew ∈Rm×1,根據本文的過程分解算法,可以將其劃分為G個子塊.每個子塊對應的時滯拓展向量可表示為xnew,g(lg)∈Rmg(lg+1)×1.需要注意的是,為方便后續描述,本文采用xnew,g表示其自身的時滯拓展向量.相應地,xnew,g所對應的4 個統計量可通過下式計算:

本文采用貝葉斯推理策略[34]融合子塊的局部評價結果,獲得全局綜合評價指標.因此,即使新的非優運行狀態所影響的過程變量分布在不同的子塊,也能獲得相對較好的評價性能.具體地,依據統計量,在線樣本子塊xnew,g處于非優狀態的后驗概率可通過式(31)計算:

式中,os和nos分別表示優運行狀態和非優運行狀態.(os)和(nos)是過程處于優運行狀態和非優運行狀態的先驗概率,可分別賦值為1-α和α.一般地,α取值為0.01 或0.05[34].

根據貝葉斯推理融合每個子塊的評價結果,全局統計量可計算為:

綜上所述,本文基于DDSFA 的運行狀態評價流程圖如圖1 所示.

圖1 基于DDSFA 的運行狀態評價流程圖Fig.1 Flow diagram of DDSFA-based operating performance assessment

3 仿真驗證

本節通過將本文的基于DDSFA 的運行狀態評價方法應用于數值仿真算例和金濕法冶金過程,以驗證其在工業過程運行狀態評價領域的有效性.為驗證分布式模型感知過程微小變化的靈敏性,將DDSFA 與集中式模型DPCA[11]和DSFA[12]進行比較.此外,為體現DDSFA 的優越性,將其與分布式模型DDPCA (Dynamic decentralized principal component analysis)[22]進行對比.

3.1 數值仿真算例

本文數值仿真算例是在文獻[9] 和文獻[22]基礎上修改得到的.數值仿真算例中包含2 個獨立的子過程.需要注意的是,雖然在實際工業過程中,較少存在像數值仿真算例這樣能理想劃分的系統,但該算例能直觀地驗證本文方法的有效性.

式中

式中,zi、ui和yi分別表示系統的狀態變量、輸入變量和輸出變量;Ai、Bi、Ci和Di表示參數矩陣;wi和ei均為零均值單位方差的隨機變量;ei的系數矩陣是為了保證yi的信噪比為1%.選擇x(t)=[y1(t),···,y4(t),u1(t),···,u4(t)]T作為過程的運行狀態評價變量.通過上述系統生成500 個樣本作為優運行狀態的訓練集,同時,在過程中引入不同的擾動,使過程偏離優運行區域生成以下2 組包含500 個樣本的測試數據.

案例1.從第101 個樣本到第500 個樣本,在w1引入幅值為2 的階躍擾動,模擬系統的潛變量發生偏移導致的非優運行狀態.

案例2.從第101 個樣本到第500 個樣本,將C1中的-0.393 改為-0.193,模擬過程變量的動態相關關系發生變化導致的非優運行狀態.

需要注意的是,上述優運行狀態和非優運行狀態是人為預定義的,并不具備實際的物理意義.

首先,將訓練數據進行標準化處理;然后,結合DTW和K-medoids 聚類算法將數值仿真算例劃分為X1=[x1,x2,x5,x6]T和X2=[x3,x4,x7,x8]T兩個子塊.不失一般性,式(7)的h取值為5.可以看出,子塊劃分的結果與實際過程的情況完全一致.針對每個子塊,通過式(19)確定其時滯系數,其中δ=1,則l1=2、l2=2.根據時滯系數將子塊訓練集進行時滯拓展,并構建對應的DSFA 模型.主導空間所保留的慢特征個數由式(22)確定.需要注意的是,在式(22)中,q=0.3.為保證對比的公平性,DPCA、DDPCA 和DSFA 都采用相同的時滯系數l=2.顯著水平α=0.01.主成分個數根據累積方差百分比(Cumulative percentage variance,CPV)確定[35],即 C PV>0.9.由于DDPCA 包含的變量子塊數等于過程變量數,且每個子塊是由不同時滯系數的過程變量構成,因此本文不作展示.

4種算法在數值仿真算例中的漏報率如表1 所示.可以看出,DPCA 和DDPCA 能有效指示過程變量的均值和方差發生的明顯變化,然而,幾乎無法感知過程變量間動態相關關系的變化.其主要原因是DPCA 采用時滯擴展矩陣的方式,將動態信息引入到原先的靜態模型中,從而實現對過程動態信息和靜態信息的提取.但DPCA 的目標函數依然主要關注過程的靜態特性,所提取的潛變量容易被靜態信息所主導,無法保證提取出準確的動態表征.因此,DPCA 和DDPCA 感知過程動態信息變化的靈敏度相對較低.DSFA 和DDSFA 對案例2 中的非優運行狀態的評價性能明顯高于DPCA 和DDPCA,其主要原因是,在建模過程中DSFA 通過T2統計量和S2統計量分別對過程變量的穩態分布p(x) 和暫態分布p(x˙) 進行單獨描述[12],從而有效地避免了過程動態信息的變化被靜態信息所掩蓋.因此,DSFA 比DPCA 和DDPCA 能更細致地描述過程數據的動態特性.需要注意的是,案例1 比案例2 所引起的變量幅值變化量大,得益于分布式建模方法更加注重過程局部信息的變化和DSFA的優良性能,DDSFA 依然獲得了令人滿意的評價結果.與DSFA 相比,DDSFA 獲得了明顯性能提升.

表1 不同算法在數值仿真算例中的漏報率(%)Table 1 Missed alarm rates of different methods in the numerical example (%)

數值仿真算例中案例1 的運行狀態評價結果如圖2 所示.從第101 個樣本到第500 個樣本,在w1中引入幅值為2 的階躍擾動,將導致過程變量x1、x2、x5、x6的均值和方差都發生明顯的變化.由于案例1 所影響的變量相對較多且變化幅值也相對較大,因而,集中式模型DPCA 和DSFA 也獲得了較好的評價結果,如圖2(a)和圖2(c)所示.由于分布式建模方法更加注重過程局部信息的變化,與集中式模型相比,分布式模型DDPCA 和DDSFA獲得了顯著的性能提升,如圖2(b)和圖2(d)所示.

圖2 數值仿真算例中,案例1 的運行狀態評價結果Fig.2 The operating performance assessment result of case 1 in the numerical example

數值仿真算例中,案例2 的運行狀態評價結果如圖3 所示.從第101 個樣本到第500 個樣本,將過程參數-0.393 改為-0.193,改變了過程變量x1、x2、x5和x6的動態相關關系.由圖3(a)、圖3(b)可以看出,DPCA 和DDPCA 并未取得令人滿意的評價結果.然而,由圖3(c)、圖3(d)可以看出,DSFA和DDSFA 評價模型中T2統計量和S2統計量均顯著地超出了控制限.得益于本文的過程分解方法將波形相近的過程變量劃分到相同的子塊,提高了DDSFA 算法對過程動/靜態信息的表征能力.因而,DDSFA 獲得了最高的評價性能.

圖3 數值仿真算例中,案例2 的運行狀態評價結果Fig.3 The operating performance assessment result of case 2 in the numerical example

3.2 金濕法冶金過程

金濕法冶金技術是從低品位礦石中提取黃金的有效方法.從生產安全角度出發,開發的算法必須經過嚴格測試才能應用于實際的生產過程.虛擬仿真技術是一種有效降低算法前期測試成本和安全風險的措施.基于以上考慮,本文課題組根據金濕法冶金過程的運行機理,設計并開發了金濕法冶金過程半實物仿真平臺.該平臺為金濕法冶金過程的監測、運行狀態評價和優化等系統的開發與驗證提供軟硬件基礎[36-37].一個典型的金濕法冶金過程包括一次氰化浸出、一次洗滌、二次氰化浸出、二次洗滌和置換5 個核心工序,金濕法冶金過程工藝流程如圖4 所示.浸出過程包含4 個自然溢流的浸出槽.向浸出槽中添加 NaCN,并充入空氣,使礦漿中的固態金溶解為 [Au(CN)2]-,此過程稱為一次氰化浸出.通過全自動立式壓濾機,實現貴液和礦漿的分離.為了使附著在濾餅中的[Au(CN)2]-盡可能少,需要使用貧液反復沖洗濾餅,此過程稱為一次洗滌.二次氰化浸出、二次洗滌與一次氰化浸出、一次洗滌類似,其作用是盡可能多地提取礦石中含有的黃金.分離出的貴液經過脫氧塔處理后,流入框板式壓濾機,并通過加入鋅粉置換出固態金.金濕法冶金過程的變量如表2 所示.

表2 金濕法冶金過程的變量Table 2 The variables of gold hydrometallurgy process

圖4 金濕法冶金過程工藝流程圖Fig.4 The flow chart of gold hydrometallurgy process

基于金濕法冶金半實物仿真平臺產生500 個優運行狀態的樣本組成訓練集.具體地,金濕法冶金過程單位時間內的綜合經濟指標位于[5 500,6 000](元/小時).同時,在過程中引入不同的擾動,使過程偏離優運行區域生成2 組包含500 個樣本的測試數據.

案例3.從第101 個樣本到第500 個樣本,將礦漿濃度由0.25 改為0.26.

案例4.從第101 個樣本到第500 個樣本,將鋅粉添加速度由0.31 線性地增加為0.33.

通過基于DTW 距離的K-medoids 聚類算法將金濕法冶金過程劃分為5 個子塊,金濕法冶金過程變量的子塊劃分結果如表3 所示,其中式(7)的h取值為5.

表3 金濕法冶金過程變量的子塊劃分結果Table 3 Sub-block division result of process variables of gold hydrometallurgy

通過式(19) 確定子塊的時滯系數,其中δ=1,則l1=l4=l5=2,l2=l3=1.顯然,本文的時滯系數確定方法根據變量子塊的特性選擇了不同的時滯系數,不僅降低了建模數據的維度,還增加了子模型的多樣性.根據時滯系數將子塊訓練集進行時滯拓展,并構建對應的DSFA 模型.主導子空間所保留的慢特征個數由式(22)確定,其中q=0.3.為了保證對比的公平性,DPCA、DDPCA 和DSFA都采用相同的時滯系數l=2.顯著水平α=0.01.主成分個數通過 C PV>0.9 確定.

4種算法在金濕法冶金過程中的漏報率如表4所示.可以看出,本文算法獲得了最高的評價性能.與DPCA 相比,DDPCA 獲得了較大性能提升.然而DDPCA 的變量子塊數等于過程變量數,當過程變量數較大時,會導致模型中的冗余信息較大,從而降低了模型性能.得益于分布式建模方法更加注重過程局部信息的變化和本文的過程分解方法提高了DDSFA 對過程動/靜態信息的表征能力,因而與DSFA 相比,DDSFA 性能獲得了明顯提升.

表4 不同算法在金濕法冶金過程中的漏報率(%)Table 4 Missed alarm rates of different methods in gold hydrometallurgy process (%)

金濕法冶金過程中,案例3 的運行狀態評價結果如圖5 所示.礦漿濃度的增加短期內會導致浸出槽中 CN-濃度下降.在 CN-濃度閉環控制器的作用下,會加大 NaCN 的添加流量.上述調節過程會導致溶液中 [Au(CN)2]-濃度升高,從而繼續影響后續工序.綜上所述,礦漿濃度變化會導致金濕法冶金過程較多的變量發生變化.因而,案例3 屬于全局非優運行狀態.由于礦漿濃度的變化量相對較小,集中式模型DPCA 并未取得令人滿意的評價結果,如圖5(a)所示.DDSFA 采用分布式建模框架,能靈敏地察覺到每一子塊的數據變化,并通過貝葉斯推理算法融合每個子塊的評價結果,因此獲得了最高評價性能,如圖5(d)所示.

圖5 金濕法冶金過程中,案例3 的運行狀態評價結果Fig.5 The operating performance assessment result of case 3 in gold hydrometallurgy process

金濕法冶金過程中案例4 的運行狀態評價結果如圖6 所示.由金濕法冶金過程工藝流程圖可知,鋅粉添加速度的變化只會影響置換過程.換言之,案例4 屬于局部非優運行狀態.由于鋅粉添加速度是線性增加的,可以清楚地看到統計量的變化趨勢.由圖6(d) 可以看出,從第270 個樣本開始,DDSFA 的統計量顯著地超出了控制限,而DSFA的統計量從第330 個樣本才開始顯著地超出控制限.綜上所述,本文DDSFA 算法對金濕法冶金過程的局部非優運行狀態具有更強的感知能力.

圖6 金濕法冶金過程中,案例4 的運行狀態評價結果Fig.6 The operating performance assessment result of case 4 in gold hydrometallurgy process

4 結束語

為充分挖掘大規模工業過程的動態信息和局部信息,本文提出一種基于慢特征分析的分布式動態工業過程運行狀態評價方法.首先,結合DTW 和K-medoids 聚類算法實現過程變量子塊的劃分,一是減少了對過程機理知識的依賴,二是充分利用了過程數據的動/靜態信息;然后,在建立子模型階段,與傳統動態建模方法對所有過程變量選擇單一的時滯系數不同的是,本文設計一種新穎的時滯系數確定方法,在降低模型的計算復雜度的同時,增加了子模型的多樣性;最后,利用貝葉斯推理融合子塊的評價結果,簡化了觸發非優運行狀態報警的最終決策.將本文方法應用于數值仿真案例和金濕法冶金過程,仿真對比結果表明,該方法具有一定的優越性.然而,本文方法主要面向線性動態工業過程,如何將其擴展到非線性動態過程或非平穩過程,值得進一步研究.