0 前言
為了提高系統的效率性、可靠性和安全性，一直以來，關於故障檢測和分離的研究備受關注。作為過程工業的一部分,化工生產過程是經過化學反應將原料轉變成產品的工藝過程。Tennessee-Eastman(TE) 過程是美國 Eastman 化學公司的 Downs & Vogel 於1993年建立的一個實際化工過程的模型,為研究該類過程的故障診斷技術提供了一個實驗平台[1]。在生產過程中,一旦發生故障而又不能及時排除時,就會影響生產的進行,並威脅到人員和設備的安全,所以,需要准確及時地識別故障類型並進行相應的處理。目前,化工過程的故障診斷方法主要有模型驅動法和數據驅動法。由於化工過程往往具有慢時變、分布參數、非線性和強耦合特性,難以得到精確的數學模型,使得模型驅動方法的性能不盡理，因此,數據驅動的方法成為處理此類問題的一個重要方向。
在工業過程一類的復雜系統中，系統規模大且過程變量高度耦合，數學機理建模難，這為故障診斷帶來了挑戰。因而，基於數學機理建模的故障診斷方法不適用於復雜系統。在當前信息時代，系統中的過程數據很豐富且容易獲取，從而數據驅動的方法被大量研究以及應用在復雜系統。其中，作為一種數據降維技術，主成分分析法（principal component analysis，PCA）能夠提取過程數據中蘊含的強相關關系。PCA 廣泛應用於復雜系統的故障診斷。如果檢測的數據點大於 PCA 模型所定義的統計限，則判斷該數據點為故障數據點，發出系統存在故障的報警。通常選取監控統計量 \(T^2\) 和 \(SPE\) 來進行故障診斷。這兩個統計量分別是基於主成分得分和殘差定義的。當 PCA 模型檢測出故障時，迸一步實施故障分離，確定哪一個或哪幾個變量是該故障的主要源頭，稱之為故障變量，從而幫助操作者根除故障源。通過假設故障變量對監控統計量具有最大的貢獻，貢獻圖（contributionplots）常用於故障分離。本文以 TE 過程作為實驗對象，應用基於 PCA 的故障診斷方法進行實驗。
1 理論
本文接下來首先介紹 TE 過程的大致工業流程和數據結構，然后在介紹 PCA 的實現原理的基礎上，給出基於 PCA 方法的故障診斷方法的實現方法步驟。
1.1 TE過程
Tennessee-Eastman(TE) 過程是 Downs 等人[1]基於 Tennessee Eastman 化學公司某實際化工生產過程提出的一個仿真系統，其具體流程見下圖。在過程系統工程領域的研究中，TE 過程是一個常用的標准問題（benchmark problem），其較好地模擬了實際復雜工業過程系統的許多典型特征，因此被作為仿真例子廣泛應用於控制、優化、過程監控與故障診斷的研究中。
TE 過程包含四種氣體原料 A、C、D 和 E，兩種液態產物 G 和 H，還包含副產物 F 和惰性氣體 B，其中進行的不可逆放熱化學反應如下：
整個過程主要包含五個操作單元：反應器、冷凝器、循環壓縮機、分離器和汽提塔。氣態反應物進入反應器，生成液態產物，反應速率服從反應動力學中的 Arrhenius 函數。產物和殘余反應物經過冷凝器冷卻后進入汽液分離器，分離得到的氣體通過壓縮機進入循環管道，與新鮮進料混合送入反應器循環使用，分離得到的液體經過管道 10 進入汽提塔進行精制，從汽提塔底得到的流股中主要包含 TE 過程的產物 G 和 H，送至下游過程。TE 過程總共包括12 個操縱變量和 41 個測量變量(包含 22 個連續過程變量和 19 個成分變量)，采用 Lyman 和 Georgakis [2] 提出的過程結構。
TE 過程模擬實際工業過程中常見的 21 種故障模式，這 21 種故障模式的具體描述如表 1 所示。本次實驗采用的TE數據集由訓練集和測試集組成，數據由 22 次不同的仿真運行數據構成，每個樣本都有 53 個觀測變量（取 52 個變量時，第 12 個操縱變量 XMV(12) 攪拌速度就不包括在內）。22 次仿真數據每次都分為了訓練集樣本和測試集樣本。正常工況下的訓練樣本是在 25h 運行仿真下獲得的，觀測數據總數為 500，而測試樣本是在 48h 運行仿真下獲得的，觀測數據總數為 960 。其余 21 次仿真數據為帶有故障的訓練集樣本，同樣分為訓練樣本和測試樣本，每組訓練樣本\測試樣本代表一種故障。帶有故障的測試集樣本在 48h 運行仿真下獲得，其中前 160 個觀測值為正常數據，故障在 8h 的時候引入，共采集 960 個觀測值。
| 故障號 | 故障描述 | 類型 |
|---|---|---|
| 1 | A／C進料比值變化，B含量不變(管道4) | 階躍 |
| 2 | B含量變化，A／C進料比值不變(管道4) | 階躍 |
| 3 | D進料溫度變化(管道2) | 階躍 |
| 4 | 反應器冷卻水入口溫度變化 | 階躍 |
| 5 | 冷凝器冷卻水入口溫度變化 | 階躍 |
| 6 | A供給量損失(管道1) | 階躍 |
| 7 | c供給管壓力頭損失(管道4) | 階躍 |
| 8 | A，B，C供給量變化(管道4) | 隨機 |
| 9 | D進料溫度變化(管道2) | 隨機 |
| 10 | C進料溫度變化(管道2) | 隨機 |
| 11 | 反應器冷卻水入口溫度變化 | 隨機 |
| 12 | 冷凝器冷卻水入口溫度 | 隨機 |
| 13 | 反應動力學參數 | 緩慢漂移 |
| 14 | 反應器冷卻水閥 | 閥粘滯 |
| 15 | 冷凝器冷卻水閥 | 閥粘滯 |
| 16 | D進料溫度變化(管道2) | 未知 |
| 17 | 未知 | 未知 |
| 18 | 未知 | 未知 |
| 19 | 未知 | 未知 |
| 20 | 未知 | 未知 |
| 21 | 閥門位置(管道4) | 卡閥 |
1.2 主成分分析
主成分分析（PCA）法也叫主元分析法，是一種數據降維算法，降維就是一種對高維度特征數據預處理方法。降維是將高維度的數據保留下最重要的一些特征，去除噪聲和不重要的特征，從而實現提升數據處理速度的目的。在實際的生產和應用中，降維在一定的信息損失范圍內，可以為我們節省大量的時間和成本。降維也成為應用非常廣泛的數據預處理方法。降維具有的優點有：使得數據集更易使用、降低算法的計算開銷、去除噪聲、使得結果容易理解。而 PCA 是一種使用最廣泛的數據降維算法。
PCA 的主要思想是將 \(n\) 維特征映射到 \(k\) 維空間中，這里的必有 \(k<n\)，在原有 \(n\) 維特征的基礎上重新構造出來的全新的 \(k\) 維特征數據具有正交特性，也被稱為主成分或主元。PCA 的工作就是從原始的高維空間中，根據方差大小順序地找出一組相互正交的坐標軸，所以新的坐標軸的選擇與原始數據密切相關。第一個新坐標軸選擇是原始數據中方差最大的方向，第二個新坐標軸選取與第一個坐標軸正交的超平面中使得方差最大的方向，第三個軸是與第 1、2 個軸正交的超平面中方差最大的，依次類推，可以得到n個這樣的坐標軸。獲得的新坐標系中我們可以發現，絕大部分方差較大的方向都包含在前面 \(k\) 個坐標軸中，其后的坐標軸所含對應的方差權重遠小於前 \(k\) 維坐標，於是只保留前面 \(k\) 個含有絕大部分方差的坐標軸。這相當於只保留包含絕大部分方差的維度特征，而忽略包含方差幾乎為 0 的特征維度，實現對數據特征的降維處理。
從以上描述我們知道得到這些包含最大差異性的主成分方向的方法，就是通過計算數據矩陣的協方差矩陣，然后得到協方差矩陣的特征值特征向量，選擇特征值最大(即方差最大)的 \(k\) 個特征所對應的特征向量組成的矩陣。這樣就可以將數據矩陣轉換到新的空間當中，實現數據特征的降維。由於得到協方差矩陣的特征值特征向量有兩種方法：特征值分解協方差矩陣、奇異值分解協方差矩陣，所以 PCA 算法有兩種實現方法：基於特征值分解協方差矩陣實現和基於 SVD 分解協方差矩陣實現。下面從最大化方差的角度簡單推導 PCA 的實現原理，並給出基於特征值分解的算法步驟。
假設有一組數據包含 \(m\) 個樣本，每個樣本 \(\boldsymbol{x_i}\in \mathbb{R}^n\)，
現在假設使用一個 \(\boldsymbol{P}\) 矩陣對樣本進行降維，
變換得到降維后的數據矩陣 \(\boldsymbol{Y}\)，其中 \(\boldsymbol{y_i}\in\mathbb{R}^k, k < n\)，
所以變換矩陣 \(\boldsymbol{P}\) 就是，
每一個 \(\boldsymbol{p}_i\) 同 \(\boldsymbol{x_i}\) 一樣是 \(n\) 維向量，那么得到這些向量就得到了 \(\boldsymbol{P}\) 矩陣。且因為由於是每一個向量與樣本進行點積的總和，所以我們應該是出於同樣的理由來確定每一個 \(\boldsymbol{p}_i\)。我們可以將每一個 \(\boldsymbol{p}_i\) 單位化，即使得 \(\boldsymbol{p^T}\boldsymbol{p}=1\)，就可以把內積 \(\boldsymbol{p}\cdot \boldsymbol{x}\) 看做是 \(\boldsymbol{x}\) 在 \(\boldsymbol{p}\) 方向上的投影，
那么為了投影過后的樣本點在投影方向上最大可能地保留原本的信息，就要讓投影點盡可能地散開，而散開的程度可以用方差來描述，所以只需要求得使方差最大的一組方向，那么矩陣 \(\boldsymbol{P}\) 中的向量就是在這些方向上的單位向量。
首先求原始樣本的協方差矩陣，
其中 \(\boldsymbol{\bar{x}}\) 是樣本均值，
則樣本點在 \(\boldsymbol{p}\) 向量所在方向上的投影的方差 \(\sigma^2\) 可以由式 \(\eqref{eq:variance}\) 表示，
其中，\(\mu = 1/m\sum_{i=1}^m{\boldsymbol{p^T}\boldsymbol{x_i}}\)，經過如下的推導過程，
由此得到如式 \(\eqref{eq:variance2}\) 的關系，
現在我們要求的 \(\boldsymbol{p}\) 向量需要滿足 \(\boldsymbol{\hat{p}}=\arg\underset{\boldsymbol{p}}{\max}{\boldsymbol{p^T \Sigma p}}\)，同時滿足前面所提的約束條件 \(\boldsymbol{p^T p}=1\)，那么我們可以通過Lagrange乘數法構造一個目標函數如式(10)，
讓目標函數對 \(\boldsymbol{p}\) 求偏導，並令偏導為零，得到，
由此可以看出，投影后的方差就是協方差矩陣的特征值。\(\Sigma\) 是一個對稱矩陣，可以做如下特征分解，
其中 \(\boldsymbol{C}\) 是一個 \(n \times n\) 的矩陣，並且是一個正交矩陣，即 \(\boldsymbol{C^T C}=\boldsymbol{E}\)。它就是由 \(\Sigma\) 的所有特征向量組成的，
前面說過，為了讓投影后的方差最大，我們需要選取 \(\Sigma\) 前 \(k\) 個最大的特征值對應的特征向量組成變換矩陣 \(\boldsymbol{P}=\begin{pmatrix}\boldsymbol{p_1}\,\boldsymbol{p_2}\cdots\boldsymbol{p_k}\end{pmatrix},\ k<n\)，從而得到降維后的數據集 \(\boldsymbol{Y}\)，滿足的關系如下
用三維變量的觀測數據舉例，如下圖所示，在三維空間中來看，這一變換的過程就是將原坐標空間變換到以特征向量方向（圖中三個箭頭所指的方向）為坐標軸的新空間中。
主元分析中需要注意的兩點是：數據的標准化和主元個數的選取問題。對於數據的標准化問題，原始數據矩陣 \(X_{m\times n}\) 每一列對應於一個測量維度，每一行對應一個樣本。兩個樣本間的相似度應正比於兩個樣本點在 \(n\) 維空間中的接近程度，由於 \(n\) 個觀測變量的量綱和變化幅度不同，其絕對值大小可能相差許多倍。為了消除量綱和變化幅度不同帶來的影響，原始建模數據應作標准化處理，
其中，\(\boldsymbol{\bar{x}_j}=1/m\sum_{i=1}^m{\boldsymbol{x}_{ij}}\) 是第 \(j\) 維上的數據均值，\(\boldsymbol{s_j}=\sqrt{1/(m-1)\sum_{i=1}^m{\left(\boldsymbol{x}_{ij}-\boldsymbol{\bar{x}_j}\right)^2}}\) 為第 \(j\) 維上的數據標准差，測試數據也要按照原始變量的均值和標准差進行標准化處理。
對於主元個數選取問題，應考慮兩個方面的因素：即原始觀測數據維數的降低和原始觀測數據信息的丟失。主元個數的選取直接影響到故障監測與診斷的效果。如果主元數目選得過小，則殘差子空間所包含的方差太多，使的殘差子空間統計量的閾值偏大，從而導致小故障難於被檢測出；而若主元數目取的太大，又會使殘差子空間包含的信息太少，使得故障對殘差影響不大，故障難於被監測出。可見，主元個數的選取是很重要的。通常，采用簡單可行的累積方差百分比（cumulative percent variance，CPV）方法來確定選取的主元的個數 \(k\)，前 \(k\) 個特征值之和除以所有特征值之和被稱為前 \(k\) 個主元的累積方差貢獻率，表示這 \(k\) 個主元所能夠解釋的數據方差的比例，即
通常選取 \(\alpha\) 為85％、90％ 或 95％. 至此我們可以得到基於特征值分解的 PCA 算法步驟：
輸入：數據集 \(\boldsymbol{X}_{m\times n}=(\boldsymbol{x_1},\boldsymbol{x_2}\cdots\boldsymbol{x_m})^T\)，需要降到k維。
- 標准化，即每一維特征減去各自的平均值，再除以每一維的標准差；
- 計算協方差矩陣 \(\Sigma_{n\times n}=\boldsymbol{X^T X}/(m-1)\)；
- 用特征值分解方法求協方差矩陣的特征值與特征向量；
- 對特征值從大到小排序，選擇其中最大的 \(k\) 個，將其對應的 \(k\) 個特征向量分別組成降維矩陣 \(P\)；
- 將數據轉換到 \(k\) 個特征向量構建的低維子空間中，即 \(\boldsymbol{Y^T} = \boldsymbol{P^T X^T}\).
1.3 基於PCA的故障檢測
為了方便描述，以下首先引入一些多元統計學中的定義。承接前面的條件，不妨假定維數 \(n=k+r\)，\(k\) 表示最終保留的維數，\(r\) 則為剩下的維數，則可以進行如下變換：
其中 \(\boldsymbol{p_i}\in\mathbb{R}^n\) 為 \(\boldsymbol{x}\) 的協方差陣特征值 \(\lambda_i\) 對應的特征向量，且各特征值如下排列
\(\boldsymbol{p_i}\) 之間是相互正交的單位向量，\(t_i\) 稱為第 \(i\) 主元，我們所關心的變量 \(\boldsymbol{x}\) 是一個 \(n\) 維向量，它的一組測量值 \(\boldsymbol{x}(\tau)\) 經過式 \(\eqref{eq:score}\) 變換后得到的值 \(t_i(\tau)\) 稱為測量值在第 \(i\) 主元上的得分。主元代表的是變量，而得分是確定值，從幾何上解釋，主元可以看做是降維后的坐標軸，而得分就是原向量在坐標軸上對應的坐標值。向量 \(\boldsymbol{p_i}\) 被稱為第 \(i\) 主元的負荷向量。
PCA 就是求取主元、得分及負荷向量的過程。我們把轉換矩陣 \(\boldsymbol{P}\) 稱為負荷矩陣，它所長成的 \(k\) 維子空間被稱為主元子空間（Principal Component Subspace, PCS），又稱模型子空間或表示子空間（Represent Subspace），\(\boldsymbol{T}=\boldsymbol{X P}\) 為由 \(k\) 列得分向量 \([t_i(1),t_i(2),\cdots,t_i(m)]^T\) 組成的矩陣稱為得分矩陣，由前面推導結論可以知道，它就是所要求的 \(\boldsymbol{Y}\) 矩陣，只是從不同角度解釋而已。后 \(r\) 個主元對應的負荷向量組成的矩陣用 \(\boldsymbol{\overline{P}}\) 表示，被稱為殘差負荷矩陣，它所張成的 \(r\) 維空間被稱為殘差子空間（Residual Subspace, RS），即是說 \(\boldsymbol{C}=\left(\boldsymbol{P,\overline{P}}\right)\)。我們可對觀測數據矩陣做如下分解，
其中的 \(\boldsymbol{\widehat{X}}_{m\times n}\) 為觀測數據矩陣 \(\boldsymbol{X}\) 在 PCS 上的估計或稱分量，\(\boldsymbol{R}_{m\times n}\) 為觀測數據在殘差子空間上的分量，被稱為模型的估計殘差矩陣。由此可知，通過對正常數據陣 \(\boldsymbol{X}\) 進行 PCA，可以把我們所關心的變量 \(\boldsymbol{x}\) 的觀測數據分解為兩個子空間的數據，
故障檢測就是檢測系統中各個監測點的數據有無異常，它通常是將新的觀測數據 \(\boldsymbol{x}_{new}\) 與系統校驗模型相比較來實現的。跟據兩者之間差距的顯著性程度，判斷系統中有無故障。主元分析法將數據空間分解為主元子空間和殘差子空間，每一組觀測數據都可以投影到這兩個子空間內。因此引入 Hotelling \(T^2\) 和平方預報誤差（Squared Prediction Error, SPE）這兩個統計量來監測故障的發生。Hotelling \(T^2\) 統計量是用來衡量包含在主元模型中的信息大小，它表示標准分值平方和。它的定義如下，
系統如果正常運行，則 \(T^2\) 應滿足
其中，\(T_\alpha\) 為 \(T^2\) 的控制限，\(k\) 為保留的主元數，\(m\) 為樣本數，\(1-\alpha\) 是置信度，\(F_\alpha(k, m-k)\) 是服從第一自由度為 \(k\)，第二自由度為 \(m-k\) 的 \(F\) 分布。\(\alpha\) 滿足概率公式 \(P\left\{F(m_1,m_2)>F_\alpha(m_1,m_2)\right\}=\alpha\)，在此定義下，\(\alpha\) 取值通常較小的值，如 0.01 左右。
\(SPE\) 統計量也稱 \(Q\) 統計量，是通過分析新的觀測數據 \(\boldsymbol{x}_{new}\) 的殘差進行故障診斷，用以表明這個采樣數據在多大程度上符合主元模型，它衡量了這個數據點不能被主元模型所描述的信息量的大小。它的計算如下，
\(SPE\) 的控制限計算公式為，
其中，
且 \(c_\alpha\) 是標准正態分布的置信極限，滿足 \(P\{N(0,1)>N_{c_\alpha}(0,1)\}=c_\alpha\)，在此定義下 \(\alpha\) 取值通常也取 0.01 左右。如系統正常運行，則樣本的 \(SPE\) 值應該滿足 \(SPE<Q_\alpha\).
在基於 PCA 的故障檢測中，對於這兩種統計量會會出現四種情況：
-
\(T^2\) 和 \(SPE\) 統計量都不超過控制限；
-
\(T^2\) 正常，但 \(SPE\) 超出 \(Q_\alpha\)；
-
\(T^2\) 超出 \(T_\alpha\)，但 \(SPE\) 正常；
-
\(T^2\) 和 \(SPE\) 都超出控制限。
一般認為 2 和 4 為故障，另兩種正常。
1.4. 基於PCA的故障隔離
當監控系統利用檢測統計量檢測到過程中發生異常后，應及時完成故障的定位與隔離（Fault Isolation），又稱故障識別（Fault Identification）。這里介紹有兩種實現故障的診斷的途徑。
一種是將故障診斷看作是模式識別的過程，所以可以采用模式識別技術直接進行故障的診斷，在這里 PCA 是作為輔助工具在條件允許時實現對數據的降維。基於模式識別的故障診斷方法包括特征提取和故障識別兩個主要步驟。第一步通過濾波消除數據信號的隨機干擾影響，提取過程特征信息，包括各種故障信息；第二步通過構建有效的分類器，將新采樣數據的過程特征同故障庫中存在的過程特征進行匹配，實現故障的診斷。常采用的模式識別技術有動態時間規整（Dynamic Time Warping, DTW）、神經網絡、支撐向量機等。但基於模式識別的故障診斷方法首先要實現過程特征的提取，這需要大量的先驗知識和含有各種故障信息的歷史數據，大大限制了該方法的應用。
另一種方法是利用基於 PCA 的技術把故障先定位到過程中某個部分或某個變量，實現故障的隔離，這樣就大大減弱了對故障庫的依賴，同時還對下一步故障診斷產生直接的指導作用，條件允許時還可以直接實現故障的診斷，在這類方法中，PCA 作為故障隔離的主要工具。貢獻圖是這類方法中的一種比較基礎的方法，一旦檢測到系統發生故障，則根據變量方差貢獻初步確定與故障關聯較緊密的若干變量或子系統，從而實現故障的定位與隔離。貢獻圖給出了被監控的各個過程變量對檢測統計量(一般為 \(T^2\) 或 \(SPE\))的貢獻，它較易生成，不需要先驗過程知識。該法的理論基礎是假設高貢獻率的過程變量是故障產生的原因，但是由於一個變量故障會影響到其它變量在主元模型下的估計值，進而影響到這些變量對檢測統計量的貢獻，所以貢獻圖只能為操作員提供一定程度的指示作用。本文的實驗以貢獻圖作為故障隔離的方法，下面給出計算過程：
當有異常狀態發生時，繪制貢獻圖，找出與故障相關的系統變量：
-
檢查每個觀測向量 \(\boldsymbol{x}\) 的標准化得分 \(t_i/\lambda_i\)，並確定造成失控狀態的 \(u\) 個得分（\(u<k\))，按下式確定，
\[\frac{t_i}{\lambda_i} < \frac{T_\alpha}{k} \] -
計算觀測向量 \(\boldsymbol{x}\) 的每個分量值 \(x_j\) 相對於失控得分 \(t_{s_i}\) 的貢獻率，按下式計算，注意下標 \(s_i\) 有 \(u\) 個而不是 \(k\) 個，它是第 1 步中確定出的 \(u\) 個得分的下標，而不是從 1 順序遞增的值，
\[cont_{s_ij} = \frac{t_{s_i}}{\lambda_{s_i}}p_{s_ij}x_j,\ i\in[1,u],s_i\in[i,k],j\in[1,n] \] -
當 \(cont_{s_ij}\)是負時，設為零，
-
分別計算每個過程分量的總貢獻率，第 \(j\) 個過程分量 \(x_j\) 的總貢獻率為
\[cont_j = \sum_{i=1}^u {cont_{s_ij}} \] -
把所有 \(n\) 個過程分量的總貢獻率畫在一個直方圖中。
2 實驗
本次實驗所用到的數據集只是全部 TE 數據集的一部分。TE過程數據集中的 22 個仿真數據集，包括 1 個正常數據集，21 個故障數據集，每個故障數據集對應一個故障模式，每個數據集都包括了 960 個采樣點，過程變量采樣間隔為 3 分鍾，成分變量的采樣間隔為 6-15 分鍾。在這 21 種故障模式中，過程起初都運行在正常模式下，在第 60 個采樣點時，發生了故障，故障一旦發生並持續到了過程結束。本次實驗采用了 52 個觀測變量，使用正常數據集中的訓練數據建立 PCA 模型，計算 \(T_\alpha\) 和 \(Q_\alpha\) 控制限，使用 21 個故障數據集中的測試集進行故障測試與隔離實驗。
2.1 實驗過程
基於 PCA 方法的在TE過程案例的故障診斷實驗流程如下圖所示。首先輸入正常數據集中的訓練集，數據規模為 \(500\times20\)。對數據標准化之后，進行 PCA 模型的建立，在建立 PCA 模型過程中，采用 90％ 的 CPV 方法選取主元數目，在得到負荷矩陣 \(\boldsymbol{P}\) 以及 \(T^2\) 和 \(SPE\) 兩種統計量的控制限之后，對故障數據集中測試集分別進行故障檢測與故障隔離實驗。
每種故障測試集的數據規模是 \(960\times52\)，在進行故障檢測時，使用建立 PCA 模型中得到負荷矩陣求取主元，而實際上得分矩陣和殘差矩陣的計算和兩種統計量的計算在可以直接通過式 \(\eqref{eq:T2}\) 和式 \(\eqref{eq:SPE}\) 在同一步驟進行計算，通過統計量與控制限的比較，判斷故障是否發生，當發生故障時，通過求貢獻圖的方式發現故障發生的可能原因，以實現故障的定位與隔離。實驗代碼（matlab）見文末。
2.2 結果分析
基於 PCA 的21種故障診斷結果如下表所示，表中粗體顯示每種故障模式下兩種統計量評判的最低故障漏報率、最低故障誤報率與最少延遲診斷時間。由於在故障模式 3、9 和 15 中過程變量數據的均值和方差均受故障影響不大，所以這三種故障難以被檢測出。從表中還可以看出，\(T^2\) 的總體誤報率要低於 \(SPE\) 統計量，而對於漏報率和延遲診斷時間上 \(SPE\) 相對有較大的優勢。所以在故障檢測過程中，可以綜合兩者的評判結果做出判斷，以得到更加准確的檢測結果。單獨從延遲診斷時間來看，基於 PCA 模型的故障診斷方法對階躍型的故障信號具有較高的敏感性，而對其他如隨機型的、閥粘滯型的或者未知類型的故障信號敏感性較弱。
| 故障號 | \(T^2\) 誤報率 / 漏報率(%) / 延遲診斷時間(T) | \(SPE\) 誤報率(%) / 漏報率(%) / 延遲診斷時間(T) |
|---|---|---|
| 1 | 0.00 / 0.63 / 4 | 8.75 / 0.13 / 1 |
| 2 | 1.25 / 1.63 / 12 | 11.88 / 0.63 / 1 |
| 3 | 1.25 / 96.88 / 20 | 18.75 / 81.25 / 1 |
| 4 | 1.88 / 45.88 / 0 | 11.25 / 0.00 / 0 |
| 5 | 1.88 / 72.63 / 0 | 11.25 / 56.50 / 0 |
| 6 | 0.00 / 0.75 / 6 | 7.50 / 0.00 / 0 |
| 7 | 0.00 / 0.00 / 0 | 9.38 / 0.00 / 0 |
| 8 | 0.63 / 2.50 / 15 | 8.13 / 2.50 / 13 |
| 9 | 6.25 / 96.38 / 0 | 16.88 / 84.88 / 2 |
| 10 | 0.63 / 54.50 / 18 | 14.38 / 27.88 / 0 |
| 11 | 0.63 / 44.50 / 5 | 16.25 / 26.50 / 6 |
| 12 | 0.63 / 1.25 / 2 | 15.00 / 1.88 / 3 |
| 13 | 0.00 / 4.63 / 36 | 8.75 / 4.00 / 17 |
| 14 | 0.63 / 0.00 / 0 | 18.13 / 1.13 / 1 |
| 15 | 0.00 / 94.50 / 91 | 9.38 / 78.88 / 2 |
| 16 | 7.50 / 74.63 / 1 | 11.88 / 34.13 / 4 |
| 17 | 0.63 / 18.38 / 1 | 20.00 / 2.50 / 10 |
| 18 | 1.25 / 10.63 / 17 | 12.50 / 7.88 / 9 |
| 19 | 0.63 / 89.25 / 10 | 8.75 / 52.88 / 1 |
| 20 | 1.25 / 59.38 / 86 | 9.38 / 30.38 / 6 |
| 21 | 3.13 / 61.13 / 26 | 24.38 / 34.63 / 0 |
從表中還可以看出基於 PCA 模型的故障檢測對TE過程的某些故障檢測效果更佳(如對故障模式1、2、4、6、12、14、17、18)，而對另外的故障情況的檢測結果不理想，這也說明了基於PCA模型的故障檢測具有局限性。
從檢測波形和貢獻圖上來看可以得到與前述結論相似的結果。基於 PCA 模型的故障診斷對於某一些故障(如下圖中的故障 14 )具有很好的檢測和隔離效果。
而對於某些故障模式，該方法有表現出不可靠的診斷效果，如下圖所示中的對故障模式 9 的診斷結果。這也說明了該方法對隨機故障信號的故障診斷效果不理想。
而從那些診斷效果較好的結果中可以看出，使用貢獻圖進行故障隔離時，使用 \(SPE\) 統計量計算的貢獻率結果相較於 \(T^2\) 統計量的計算結果更具有參考意義，如下圖中第一張圖中所示，其計算出的 52 觀測變量對於當前故障的貢獻率大小區分度更為明顯。
而基於貢獻圖的方法也被證明了具有很大局限性，Alcala 和 Qin[4] 證明貢獻圖對某些簡單的傳感器故障也不能保證得到正確的故障診斷結果，並為此提出了基於重構的貢獻圖方法(reconstruction based contribution，RBC)。更為本質的原因是像 PCA 這類的純數據驅動方法無法描述過程的機理反應，所建的數據模型缺乏清晰的物理含義，相互之間的因果關系不明確，因此當應用在復雜系統中定位和診斷故障源時就造成了很大限制。同時基本的 PCA 模型沒有考慮到很多更復雜的問題，如非線性問題、局部性問題、稀疏性問題、動態性問題、自適應問題等，這些實際應用中的問題也是其診斷效果十分局限的原因[7]。
實驗代碼
% 基於 PCA 模型的在 TE 過程中的故障診斷仿真實驗
close all; clear; clc;
data=load('TE_data.mat');
data = struct2cell(data);
testdata = data(1:22);
train = cell2mat(data(23));
train = train';
train_mean = mean(train); %按列 Xtrain 平均值
train_std = std(train); %求標准差
[train_row,train_col] = size(train); %求 train 行、列數
train=(train-repmat(train_mean,train_row,1))./repmat(train_std,train_row,1); %標准化
sigmatrain = cov(train); %求協方差矩陣
%對協方差矩陣進行特征分解，lamda 為特征值構成的對角陣，T的列為單位特征向量，且與 lamda 中的特征值一一對應：
[T,lamda] = eig(sigmatrain);
%取對角元素(結果為一列向量)，即 lamda 值，並上下反轉使其從大到小排列，主元個數初值為 1，若累計貢獻率小於90%，則增加主元個數
D = flipud(diag(lamda));
num_pc = 1;
while sum(D(1:num_pc))/sum(D) < 0.9
num_pc = num_pc +1;
end
P = T(:,train_col-num_pc+1:train_col); %取與 lamda 相對應的特征向量
%每一列代表一個特征向量
%求置信度為 99%、95%時的 T2 統計控制限
T2UCL1=num_pc*(train_row-1)*(train_row+1)*finv(0.99,num_pc,train_row - num_pc)/(train_row*(train_row - num_pc));
T2UCL2=num_pc*(train_row-1)*(train_row+1)*finv(0.95,num_pc,train_row - num_pc)/(train_row*(train_row - num_pc));
%開始計算SPE統計量
theta = zeros(3,1);
for i = 1:3
theta(i) = sum((D(num_pc+1:train_col)).^i);
end
h0 = 1 - 2*theta(1)*theta(3)/(3*theta(2)^2);
ca = norminv(0.99,0,1);
SPE = theta(1)*(h0*ca*sqrt(2*theta(2))/theta(1) + 1 + theta(2)*h0*(h0 - 1)/theta(1)^2)^(1/h0);
frecord = fopen('./records.txt','w');
for k = 1:21
%第k組測試數據
fprintf(frecord,'故障模式(%d)—',k);
test = cell2mat(testdata(k+1));
%開始在線檢測
n = size(test,1);
test=(test-repmat(train_mean,n,1))./repmat(train_std,n,1);
%測試樣本歸一化
[r,y] = size(P*P');
I = eye(r,y); %單位矩陣
T2_test = zeros(n,1);
SPE_test = zeros(n,1);
T2_falm = 0; SPE_falm = 0;
xf_T2 = []; cf_SPE = [];
for i = 1:n
T2_test(i)=test(i,:)*P/lamda(train_col-num_pc+1:train_col,train_col-num_pc+1:train_col)*P'*test(i,:)';
SPE_test(i) = test(i,:)*(I - P*P')*test(i,:)';
if T2_test(i) > T2UCL1
if i < 161
T2_falm = T2_falm + 1;
else
xf_T2 = cat(1, xf_T2, i);
end
end
if SPE_test(i) > SPE
if i < 161
SPE_falm = SPE_falm + 1;
else
cf_SPE = cat(1, cf_SPE, i);
end
end
end
fprintf(frecord,'T2,SPE(誤報數/誤報率/漏報率/延遲診斷時間)：%3d/%.2f/%.2f/%d, %3d/%.2f/%.2f/%d\n',T2_falm,T2_falm/1.6,100-length(xf_T2)/8,xf_T2(1)-161, SPE_falm,SPE_falm/1.6,100-length(cf_SPE)/8,cf_SPE(1)-161);
%繪圖
figure;
% suptitle(['故障模式',num2str(k-1)]);
subplot(221);
plot(1:n,T2_test,'k');
title('主元分析統計量變化圖T2');
xlabel('采樣數');
ylabel('T^2');
hold on;
line([0,n],[T2UCL1,T2UCL1],'LineStyle','--','Color','r');%畫出標志線
line([0,n],[T2UCL2,T2UCL2],'LineStyle','--','Color','g');
subplot(223);
plot(1:n,SPE_test,'k');
title('主元分析統計量變化圖SPE');
xlabel('采樣數');
ylabel('SPE');
hold on;
line([0,n],[SPE,SPE],'LineStyle','--','Color','r');
%貢獻圖
%1.確定造成失控狀態的得分
pchs = randsample(cf_SPE,1);
fprintf(frecord,'用於計算貢獻率的故障點：%d\n',pchs);
S = test(pchs,:)*P(:,1:num_pc);
r = [ ];
for i = 1:num_pc
if S(i)^2/lamda(i) > T2UCL1/num_pc
r = cat(2,r,i);
end
end
%2.計算每個變量相對於上述失控得分的貢獻
cont = zeros(length(r),52);
for i = 1:length(r)
for j = 1:52
tmp = abs(S(r(i))/D(r(i))*P(j,r(i))*test(pchs,j));
if tmp < 0
tmp = 0;
end
cont(i,j) = tmp;
end
end
%3.計算每個變量的總貢獻
CONTJ = zeros(52,1);
for j = 1:52
CONTJ(j) = sum(cont(:,j));
end
%4.計算每個變量對Q的貢獻
e = test(pchs,:)*(I - P*P');
contq = e.^2;
%5. 繪制貢獻圖
% figure;
subplot(222);
bar(CONTJ*100,'k');
title('T^2貢獻率');
xlabel('變量號');
ylabel('%');
subplot(224);
bar(contq*100,'k');
title('Q貢獻率');
xlabel('變量號');
ylabel('%');
end
fclose(frecord);
參考
-
JJ. Downs,E.F. Vogel, A plant-wide industrial process control problem[J]. Computers & Chemical Engineering, 1993. 17(3): p. 245-255.
-
P.R. Lyman,C. Georgakis, Plant-wide control of the Tennessee Eastman problem[J]. Computers & Chemical Engineering, 1995. 19(3): p. 321-331.
-
B.M. Wise,N.B. Gallagher, The process chemometrics approach to process monitoring and fault detection[J]. Journal Of Process Control, 1996. 6(6): p. 329 348.
-
C.F. Alcala, S.J. Qin, Reconstruction-based contribution for process monitoring[J]. Automatica, 2009. 45(7): p. 1593-1600.
-
Ian T. Jolliffe. Principal component analysis : A beginner’ s guide[J]. Weather, 1993. 48(8): p. 246–253.
-
Ian. T. Jollife, J. Cadima. Principal component analysis: A review and recent developments[J]. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2016. 374(2065).
-
劉康玲. 基於自適應PCA和時序邏輯的動態系統故障診斷研究[D].浙江大學,2017.