優化動態核主元分析的工業過程故障監測方法

楊 芳， 王亞君， 沈亞慧

（遼寧工業大學電子與信息工程學院，遼寧錦州 121001）

0 引 言

故障監測一直是研究熱點。復雜工業過程往往具有非線性、動態、多模態、高維度、間歇性等特點，使得故障監測變得非常復雜。傳統的故障監測方法已經無法滿足故障監測的要求，需要采用更加先進的技術實現更準確地故障監測［1-3］。Scholkopf等［4］提出了核主元分析（KPCA）方法，該方法使用核函數實現某種非線性映射，解決了之前存在的線性映射無法解決的問題。Chen等［5］提出了BDPCA（batch dynamic PCA）和BDPLS（batch dynamic PLS）算法，然而擴展數據陣的有效范圍受時滯長度的限制，僅能反映局部時間段內的動態特性。Wang等［6］提出了一種基于多動態核主元分析（MDKPCA）的復合降維故障監測方法。將自回歸移動平均時間序列模型和KPCA 方法相結合，建立動態核主元分析（DKPCA）模型，對工業過程中的非線性和動態進行處理。

目前KPCA方法作為一種非線性主元分析方法被廣泛運用，KPCA 方法的有效性依賴核函數參數的選取，因此關于核函數參數的研究一直是近幾年的熱點。Chapelle等［7］提出了利用梯度下降法解決核函數優化問題，但需要計算目標函數對優化參數的偏導數。如果目標函數對某個參數的偏導數不存在，或者由于計算復雜而難以求解，梯度下降法就不能有效地進行核函數優化。Zhang等［8］進行了核函數參數優化方法的比較研究，其中包括二分法、基于反向傳播（BP）神經網絡的重構法以及基于樣本分類的重構法，并設計了一個核函數參數優化系統進行核函數參數優化，但存在計算量較大的問題。Chen 等［9］提出了一種改進的遺傳算法，用于核函數優化。將核函數的類型和參數作為決策變量，并通過優化編碼方案、遺傳策略、適應度函數、交叉和變異算法等方法最小化均方根誤差。為了確保核函數的多樣性，采用二折交叉驗證法對訓練結果進行驗證，但該方法容易陷入局部最優。

針對上述問題，提出了一種基于人工大猩猩部隊優化動態核主元分析（GTO-DKPCA）的故障監測方法。首先將數據進行預處理，然后利用時間序列模型和KPCA方法建立故障監測模型，最大化人工大猩猩部隊優化算法對正常數據和故障數據之間的幾何距離來尋找最優參數。為了驗證該方法的有效性，給出了青霉素發酵故障監測工藝案例。

1 動態核主元分析

1.1 核主元分析過程

KPCA方法主要解決非線性問題，基本思想是：使用非線性映射將原輸入空間映射到一個髙維的特征空間F中，然后在該特征空間中進行主成分分析計算，從而把原始空間的非線性問題轉化為高維空間的線性問題。KPCA方法在提取數據非線性主要成分的同時保留了原始數據結構，很好地解決了非線性特征變換問題，也不會損失數據的重要特征。KPCA方法如下：

設樣本數據xk∈Rm（k＝1，…，N），其中N為樣本點個數。令F為特征空間，將樣本數據映射到特征空間：xk→F，映射后的數據為Φ（xk）。協方差矩陣計算式為

對應特征方程為

式中：λ為非零特征值；v為特征向量。式（2）等價于

存在系數αi（i ＝1，2，…，N），使得特征向量v＝結合式（2）、（3）得

為了避免特征空間計算復雜，引入核函數K，化簡可以得到Nλα ＝Kα，其中K滿足中心化公式INK－KIN＋INKIN，IN＝1／N。計算過程詳見文獻［10］。第m個KPCA的特征主元通過映射得到，計算過程為

式中：um為非線性主元；q為非零特征值個數；vm為對應的特征向量；αi，m為um對應的系數；m＝1，2，…，q。

1.2 動態核主元分析過程

為了解決生產過程中的動態特性，構造動態自回歸移動平均時間序列模型［11］。當前時刻每個采樣點的變化主要與時間序列的歷史數據有關，利用時間序列歷史相關特性在每個采樣點處增加前d個歷史時刻上的采樣點。時間序列模型與KPCA方法相結合稱為DKPCA方法。DKPCA方法可以很好地去除動態系統中輸入與輸出的關系，能夠更好地監測故障。擴展后的數據矩陣為：

式中：Xd為擴展數據矩陣；H為每個批次建模數據矩陣的采樣點；d為時滯長度。

1.3 相關統計量與控制量計算

DKPCA方法的統計量和控制量與PCA方法的統計量類似，但對于DKPCA方法，需要在特征空間中進行計算。在特征空間中平方預測誤差（SPE）統計量計算方法為

假設預報誤差服從χ2分布，此時特征空間中SPE統計量的控制限可由式（7）計算得到。在特征空間中，霍特林統計量T2計算方法也與PCA方法類似，具體計算方法為

式中，Λ－1為對角陣的逆矩陣。對角陣由與um相對應的主元方差構成，故特征空間中T2控制限可由式（8）計算獲得。

2 基于人工大猩猩部隊優化算法的核函數參數優化

2.1 核函數選取

將核函數的內積代替復雜的非線性映射來處理非線性問題。通過核函數可以將測量樣本數據映射到高維特征空間，在高維特征空間中應用線性主元分析，因此核函數的選取尤為重要。由文獻［12］可知，多項式核函數、徑向基核函數（RBF）和Sigmoid 核函數是常用的核函數。研究表明［13-15］，徑向基核函數對數據的處理更好，因此選擇徑向基核函數。徑向基核函數表達式為

式中：c為核函數參數，表示徑向基核函數的寬度。在進行非線性特征變換時，核函數參數的選取對DKPCA方法的監測性能有著重要的影響。合適的核函數可以對數據進行更精準的特征提取，有效地避免了非線性情況下數據不可分的難題，從而提高故障監測的準確性和可靠性。傳統的核函數參數一般由人為主觀測試或交叉驗證等方法選取。人為選取用時較長，而交叉驗證方法需要大量的計算。因此，提出了GTO-DKPCA方法，利用優化算法對核函數參數進行優化，得到最優參數。

2.2 人工大猩猩部隊優化算法

元啟發式算法在解決優化問題時起著至關重要的作用，GTO算法靈感來自大猩猩在自然界中的社會行為［16］。在GTO 算法的初始化過程中先設置相應參數，然后在探索空間中隨機初始化種群。Q為種群數量，M為最大迭代次數。β 為大猩猩暴力程度公式中的參數，其值為0.3。w是開發階段2 種機制的控制參數，其值為0.8。探索階段模型為：

式中：G（t＋1）為大猩猩下一次迭代時的位置；Y（t）為當前大猩猩的位置；r1、r2、r3和rand 為每次迭代更新范圍0 到1 的隨機值；p為介于0 到1 之間的給定參數；U、L為變量的上、下限；Yr（t）和Gr（t）為隨機選取大猩猩的位置；t為當前迭代次數。在每次更新位置分別計算D、B、R，計算式為：

式中：r4為迭代時更新從0 到1 的隨機值；l為0 到1之間的隨機數；Z為區間［－D，D］內的隨機數。在探索階段結束時，計算G（t）和Y（t）的適應度，如果G（t）適應度小于Y（t），則Y（t）被G（t）代替。因此，該階段產生的最優個體也被視為銀背大猩猩。在GTO 算法的開發階段，模擬了大猩猩的2 種行為：跟隨銀背大猩猩和成年雌性競爭。這2 種行為通過D進行切換，當D≥w時，選擇跟隨銀背大猩猩，相反則進行競爭。跟隨銀背大猩猩機制的表達式為：

式中：Ysilverback為最優解；Ga（t）為迭代時每只候選大猩猩的位置。成年雌性競爭機制的表達式為：

式中：I為大猩猩競爭強度；r5為0 和1 之間的隨機數；A為用來模擬競爭程度的系數；O為暴力對解決方案維度的影響。rand為0 到1 之間的隨機數。當rand≥0.5 時，O等于正態分布和問題維度中的隨機數；反之，O等于正態分布中的隨機數。人工大猩猩部隊優化算法模擬如圖1 所示。

圖1 人工大猩猩部隊優化算法模擬

2.3 基于GTO-DKPCA故障監測方法

采用GTO方法優化徑向基函數的參數c，使用自適應度函數進行故障監測。給定數據集，計算正常數據集與故障數據集之間的類間幾何距離與數據點到最近類距離的比值，目標是最小化距離比值的平均值。該自適應度函數能夠適用于正常數據集和故障數據集，從而提高故障監測的準確性。該適應度函數定義為

式中：Sd為數據間的散度矩陣；Sb為數據中的散度矩陣。根據正常數據與故障數據之間的關系，建立適應度函數，利用適應度函數對核函數參數尋優，更好地平衡正常數據和故障數據之間的差異，提高故障監測的準確性和可靠性。

3 實驗驗證

青霉素發酵生產過程是一個典型的非線性、動態性和間歇性過程。在仿真實驗中對通風率、攪拌功率和底物供給速率3 種常見故障類型進行研究。青霉素發酵過程中，前2 種故障主要影響溫度和溶解氧濃度，從而影響青霉素的產量和質量。正常工況下采集30批次建模數據作為訓練數據。各批次的反應時間設定為400 h，進行數據采集時采樣間隔為0.05 h，建模和在線監測時采樣間隔為1.00 h。初始條件和設定值范圍如表1 所示。

表1 發酵過程初始參數取值

在青霉素發酵過程中，30 批次的訓練數據中包含10 個變量8 000 個數據，在利用GTO算法對徑向基核函數參數進行優化時，對正常數據與故障數據的類間幾何距離最大適應度函數尋優，最終得到優化后參數。采用GTO算法進行優化的適應度曲線如圖2 所示。

圖2 適應度曲線

為了驗證所提方法的有效性，在正常工況下采集30 批次建模數據，分別建立多向核主元分析（MKPCA）、BDKPCA、GTO-DKPCA 模型。采用置信水平為99%的控制線進行在線監測。表2 列出了青霉素發酵過程的3 種測試情況：攪拌功率階躍跳變、通風率線性變化、底物供給速率線性變化。

表2 青霉素發酵過程中的3 種模擬場景

圖3 為故障1 監測結果。在故障1 中，通風率在200 h時以0.05 的斜率線性降低，直至反應結束。對于MKPCA方法（見圖3（a）），T2統計量在211 h時發生偏差，SPE統計量在207 h時發生顯著偏差，T2統計量故障監測延遲了11 h，SPE 統計量故障監測延遲了7 h。BDKPCA方法的監測結果如圖3（b）所示，SPE統計量從205 h開始監測到故障，而T2統計量從209 h開始監測到故障。GTO-DKPCA 方法能夠在203 h 時監測到故障。MKPCA和BDKPCA方法不能夠很明顯地監測出異常情況，而GTO-DKPCA 方法能夠很好地監測出異常，如圖3（c）所示。當干擾在200 h 引入時，相較于MKPCA 和BDKPCA 方法，GTO-DKPCA 方法監測能力更強。

圖3 故障1監測結果

由于底物供給速率的快慢會對青霉素產量產生影響，因此將底物供給速率在50 h 時增加10%至結束。如圖4（a）所示，MKPCA方法進行故障監測時，SPE 和T2統計量均在50 h之前就監測出故障，將正常數據作為故障數據；如圖4（b）所示，BDKPCA 方法進行故障監測時，T2統計量在62 h 監測出故障，SPE 統計量在39 h監測出故障；如圖4（c）所示，GTO-DKPCA方法監測故障時，T2統計量在59 h 時監測出故障，SPE 統計量在51 h時監測出故障。相較于MKPCA與BDKPCA方法而言，GTO-DKPCA方法能夠更好地監測故障。

圖4 故障2監測結果

對于故障3，攪拌功率從100 h 開始階躍增加10%直至反應結束。攪拌功率變化將直接影響總體氧氣質量系數K1a，而K1a的降低導致培養基中溶氧水平降低，進而降低了生物量增長和青霉素濃度。圖5（a）、（b）分別為MKPCA 和BDKPCA 方法的監測結果。MKPCA 無法有效地監測出故障，而對于BDKPCA，在101 h時SPE統計量監測出故障，但T2統計量仍不能監測出故障。GTO-DKPCA 方法能夠成功地監測出故障，如圖5（c）所示。

圖5 故障3監測結果

4 結 語

提出了一種GTO-DKPCA 故障監測方法。首先，根據正常數據和故障數據的特點，利用GTO方法對核函數參數進行優化，使DKPCA模型的監測效果更好。通過計算特征空間中的T2和SPE統計量，有效地提取過程變量之間的非線性關系，從而實現更佳的監測性能。在青霉素發酵過程中應用GTO-DKPCA方法進行故障監測。監測結果表明，GTO-DKPCA 方法比MKPCA和BDKPCA方法更準確地監測出故障，具有更為優越的監測性能。