• 沒有找到結果。

利用動態數據做系統之監控偵錯與數據調諧(2/3)

N/A
N/A
Protected

Academic year: 2021

Share "利用動態數據做系統之監控偵錯與數據調諧(2/3)"

Copied!
12
0
0

加載中.... (立即查看全文)

全文

(1)

行政院國家科學委員會專題研究計畫 期中進度報告

利用動態數據做系統之監控偵錯與數據調諧(2/3)

計畫類別: 個別型計畫

計畫編號: NSC93-2214-E-002-014-

執行期間: 93 年 08 月 01 日至 94 年 07 月 31 日

執行單位: 國立臺灣大學化學工程學系暨研究所

計畫主持人: 黃孝平

計畫參與人員: 鄭智成 林大溱

報告類型: 精簡報告

處理方式: 本計畫可公開查詢

中 華 民 國 94 年 6 月 1 日

(2)

利用動態數據做系統之監控偵錯與數據調諧 (2/3)

Abstract

The objective of this research is to utilize dynamic-PLSR1 (DPLS) to identifying dynamics of a process, and utilize PLS on modeling wavelet transformed input data for FDI. After modeling the dynamics of the process, filters are constructed and applied to the input data which have been pre-treated with Discrete Wavelet Transform (DWT). The DWT decomposes each input signal into three basic components -- trend, seasonal, and, stationary and random components. By this DWT pretreatment, it provides a capability to identify the source of sensor faults. The concepts and techniques are demonstrated using a simulated Wood and Berry system. It shows this presented method is effective in detecting and identifying either the faults caused by high frequency noises or by biases in relating outputs.

Keywords: s

ubspace identification, impulse response sequence, closed loop, PLS, fault detection,

wavelets

一、 緣由與目的

製程偵錯與診斷(Fault Detection and Isolation, FDI)的目標,為儘早發現製程操 作異常問題並了解造成異常的原因及時 間。透過製程的監測,以維持製程操作的 理想狀況,避免不必要的產品損失與成本 的浪費。亦即提高產品的品質與良率,發 揮製程的最大效能。 多 變 量 製 程 管 制 技 術 (Multivariate statistical process control, MSPC)由於能克 服 諸 多 品 質 與 製 程 變 數 共 線 性 (Collinearity)的問題[1],被大量應用於產 品品質變數監控與可量測的製程變數管 制與監測。此乃因大量多變數間若存在內 相依性(inter-dependence),則單一量測變 數的變化不能作為正確的參考指標。因 此,早期的多變量統計管制圖,利用 T2 及 Q 管制圖的統計計算原理,克服變數共 線性的問題。後期更結合多變量投影法 (Multivariate Projection Method)原理,如主 變 量 投 影 法 (Principal Components Analysis, PCA)及部分最小平方法(Partial Least Squares, PLS),將諸多可能相關之變 數構成的複雜向量空間投影至較低維度 且完全線性獨立的潛藏變數空間(Latent variable spaces),以擷取大量變數產生龐 大數據隱含的重要訊息。透過模型的建立 與驗證,對正常與異常操作狀態作更明確 的界定,在製程品質監控上佔有一席之 地。 化工製程的偵測,必須透過量測器 (Sensor) 取 得 製 程 變 數 及 輸 出 的 正 確 訊 息,才能了解製程異常及錯誤(fault)發生 的精確時間及其型態。但是工廠的量測器 常有校正異常及雜訊(noise)干擾過大的問 題,且大部分的重要訊息發生於短暫的時 間之內。因此,小波分析(Wavelet Analysis) 提供訊號分解及去除雜訊(de-noising)[16] 的工具。Kosanovich 及 Piovoso(1997)[17] 曾應用小波轉換的優點,結合前人以 PCA 及 PLS 在製程偵測的方法。 本計劃應用一般的多變量製程管制 技術(Multivariate statistical process control) 與偵錯技術,配合傳統製程識別的概念以 較為簡化的技術建立製程模型,以此模型 預測輸出端是否有錯誤發生及輸出端開 始發生錯誤的時間點。並於輸出端執行錯 誤診斷,除貢獻管制圖的使用,並配合小 波轉換(Wavelet Transform)分解訊號變化 的型態,更進一步診斷錯誤的來源,稱之 為 預 測 式 錯 誤 診 斷 (Predicted Fault diagnosis),如圖一。 Process Model (T. F. filter ) Æ Filtered Input PLS2 Model R R1 R2 R3 S S1 S2 S3 ˆD X ˆ B X D X B X Process Information Predicted Faults Diagonalsis Wavelet Decomposition PLS Monitoring 圖 一、預測式錯誤診斷(Predictive Fault Diagnosis)示意圖

(3)

二、 部分最小平方法(PLS)

近代工業製程中,由於大量複雜的 製程變數且常隱含及高的相關性。MSPC 由於可以擷取其中的重要訊息,並將維 度作大幅簡化而受到廣泛的應用。其中 部分最小平方法(Partial Least Squares, PLS)就是常見的一種重要工具之一。

部分最小平方法(Partial Least Squares, PLS),亦稱為潛結構投影法(Projection to Latent Structure) , 首 先 由 經 濟 學 家 Wood(1966)[25]所提出,初期廣泛為心理 學等社會科學領域,化學與藥學等自然科 學等領域所使用。PLS 之理論與方法,近 幾年受到數學家與各方的質疑與討論。發 展 至 1980 年 以 後 , 開 始 有 非 線 性 應 用 (Nonlinear PLS) 、 Prefiltered PLS 、 Orthogonal Signal Correction (OSC) 、 Orthogonalized-PLS (O-PLS)等的發展。並 於時間序列(time series)、批次(batch)及 wavelets variants。 PLS 與 傳 統 多 變 量 最 小 方 差 分 析 (MLR) , 皆 是 處 理 兩 群 變 數 之 相 關 性 (relation)問題。但前者的發展與應用潛力 更優於後者,主要因為 PLS 隱含潛變數 (latent variable)的概念及克服多變量所衍 伸的共線性問題。且適用於現實生活中經 常 發 生 的 非 因 果 關 係 模 型 (non-causal model)情況。PLS 引用 PCA 的概念,分析 兩 組 可 能 各 自 內 相 關 的 多 變 量 數 據 X ( predictor block )與 Y ( predicted block )之 潛變數,如(2)及(3)式。但由於必須找到兩 組數據之主要關聯性的實際需求,而將 X 與 Y 合併成為(4)式。 a T T j j j=1 a T T j j j=1 a T T j j j=1 X = TP + E = t p + E Y = UQ + F = u q + F ˆ ˆ ˆ ˆ Y = UQ + F = u q + F, (U = TB)

(2) (3) (4) 由於計算方式及物理意義的不同, PLSR1 由 於 數 學 意 義 較 為 簡 單 而 實 用 (Helland, 2001))[29]。因此,在識別及預測 的應用上較為廣泛。PLSR1(PLSR of One predicted variables),每次運算僅考慮單一 被預測變數 yi;如此,選取的變量僅與此被 預測變數有最大共相關性。對於多被預測 變量系統,以一連串的 PLSR1 演算法分別 對每一被預測變數作 PLS 分析,而得到多 組分別針對不同被預測變數的 X 變量。 1 p 0 i,0 i T j-1 i i, j T j-1 j (Y = y ... y ) ( i, j = 0 , E = X , f = y ) E y (1) w = E y     ⋅ ⋅ 起始值 i, j j-1 i, j T j-1 i, j i, j T i, j i, j

j,new j,old j,old

j,new j,old j,old

j,new j,old j,old

T j j-1 i, j i, j T i i, j i, j i, j-1 i, j i, j i, j j T i, j i, j (2) t = E w E t (3) p = t t (4) p = p p (5) t = t p (6) w = w p (7) j = j +1, to j = min(m, n) E = E - t p y t f = f - b t q , b = t t ⋅ ⋅ ⋅ ⋅ ⋅ 同理可將新數據投影至建模完成的主軸 空間(PCS),並依照 X 和 Y 的關係預測 Y 值,可估計如下: 1,a a ˆ Y = X B1⋅ (6) 其中 a 1, j 2, j p, j T -1 T -1 T i, j i, j i, j i, j i, j i, j i, j 0, j B1 = b b b , b = W (P W ) (T T ) T f , j a     = L

三、 PLS 在偵錯上的應用

另一方面,以多變量投影法(如 PCA 及 PLS 等)結合 Hotellin T2及 Q 管制圖於 統 計 製 程 偵 測 (Statistical Process Monitoring, SPM)的使用,必須符合時間獨 立(Time independency)且量測變數與 PCA 的殘差(Residuals)需為常態分布(Normal distribution)假設的假設條件。然而,實際

(4)

的化學製程中通常包含動態(dynamic)及 非線性(non-linearly)的關係。因此,wise 等 人 (1990)[18] 提 出 動 態 自 相 關 (auto-correlation)發生時,通常不滿足常態 分布的假設。於是 Ku (1995)等人提出動態 主 成 分 分 析 法 (Dynamic Principal Component Analysis, DPCA)[3]。將時延 (time lag) 概 念 的 造 成 的 動 態 關 係 加 入 PCA 的模型之中,以克服動態及自相關性 (autocorrelation)在統計處理的限制。 PLS 在偵錯上的應用上,由於其探討 兩類(two-block)變數變數:預測因子 X 及 被預測因子 Y 之關聯性,因此 PLS 在動 態或製程輸入及輸出產生相關變數的診 斷 技 術 上 較 PCA 為 多 元 (Multiple) 。 MacGregor 等(1995)[15]強調 PLS 在工程 偵測應用時,除可利用 X 變量(t-score)的 Hotelling T2管制圖外,亦可分析 X 及 Y 變 數 個 別 的 平 方 預 測 誤 差 (Squared prediction error, SPE )。Kruger (2001)[12] 更提出可以對 Y 變量(u-score)作平方預測 誤 差 分 析 , 以 檢 視 潛 變 數 模 型 (latent variable model)解釋製程變異的程度。 1. 動態 PLS 及其於識別的應用 同理,PLS 在動態模型的識別上, Ricker(1988)[8]首先提出將動態時延參數 加入輸入資料計算 PLS 模型,並可轉換成 PLS 有限脈衝應答(Finite Impulse Respond, FIR)並識別為轉移函數(Transfer Function) 模型。另外,Qin and McAvoy (1992a)[13] 提出另一將輸入及輸出都加入時延的方 法 , 可 以 得 到 ARMA(Autoregressive Moving Average)模型。另外,為克服非線 性 的 問 題 , Qin and McAvoy (1992b,1996)[34,40]提出積分 PLS 類神經 網 路 結 構 (Integrated PLS-neural net structure)。其中前兩識別方法牽涉時延階 次的選定,而使過程繁複而難以處理。因 此,Kaspar and Ray(1992)[9,10]提出另一 動 態 轉 換 (Dynamic Transformation) 的 方 法,利用動態濾波器(dynamic filter)將輸入 時間序列作轉換,使輸入和輸出間僅存線 性關係。如此,可直接對轉換過的變數作 線性代數(Linear Algebra Relation)分析,但 濾波器的階數必須審慎評估以克服模型

維度(automatic model reduction)不在簡化 的缺點。 上述動態識別的發展,主要應用於批 次(batch type)識別,亦即須先將數據收集 完成後予以建模及偵錯。如此,對於動態 模式改變或大量數據貯存所需之大量空 間之衍生之問題,有其限制及缺點。因 此 , Ljung (1911)[14] 提 出 遞 迴 PLS(recursive PLS)的方法,線上(on-line) 更新動態模式。後續,Qin(1993, 1998)對 化工製程得應用及並提出各式遞迴 PLS 的 演 算 法 。 本 計 劃 於 程 序 識 別 使 用 Ricker(1988) 方 法 並 比 照 動 態 PCA (Ku,1995)[3],稱之為動態偏最小回歸法 (dynamic PLSR)。 2 動態 PLSR1 及程序識別 由於 PLSR1 僅針對單一預測輸出 (predicted output)的理論結構,相較於 PLS2 為簡單。因此,實際演算上不需處 理多預測輸出的複雜關係,模型的精確 度較高。所以,仿效 DPCA 改寫 PLSR1 的樣本矩陣排列如下,以計算動態系統 的 脈 衝 應 答 係 數 (impulse response coefficients) , 而 可 對 動 態 系 統 作 識 別 (Identification)。 每次取單一輸出變數將樣本矩陣排 列如下: 1 1 1 (1) (1) (2) (2) ( ) ( ) u u u n n n x x x x x x N x N       =       L L M M L (1) (2) ( ) y y y y N       =     M (7) (8) 其中

x

為輸入矩陣,包含

n

u個預測變數及 N 個取樣,x∈RN×nu。而y為輸出矩陣, 包含單一個被預測變數的 N 個取樣點,

(5)

N 1 R y

×。在執行 PLSR1 之前,同樣對原 數據均值原點化及正規化的前處理,若給 定 n 則可重新排列矩陣如下: 1 1 1 1 1 1 ( ) (1) ( 1) (2) ( 1) ( ) ( ) (1) ( 1) (2) ( 1) ( ) u u u u u u reg n n n n n n x n x x n x x x N x N n x n x x n x x N x N n   +  =    +   − −  L L L L M M L L L L L L M M L L ( 1) ( 2) ( ) reg y n y n y y N +    +    =     M (9) (10) 其中 n 為脈衝應答係數的係數個數。xreg 為回歸 X 矩陣(Regression X Matrix), (N-n-1) ( )

R

n nu reg

x

× × 。yreg為回歸 Y 矩陣 (Regression Y Matrix),

y

reg

R

(N-n-1) 1×。如 此,以上述矩陣排列執行 PLSR1,稱之為 動態 PLSR1 (dynamic-PLSR1)。將所得之 回歸係數序列(亦即 3.2 述及之 PLSR1 演 算法中(6)式的 B1a或 b1,j),再以相對應的 不同預測變數分開計算。則可得到每一預 測變數相對於此被預測變數的脈衝應答 係數(以

θ

表示)和 PLSR1 回歸係數的關係 如下:

[

b1a(1: )n b1a(n 1: 2 )n b1a( : )l r

]

n m θ= + L × (11) 其中

l

= − × +

(

n

u

1)

n

1,

r n n

= ×

u 。將所得之脈 衝應答係數,經過累加的計算可以得到其 階梯應答(Step Respond)曲線。如此,可由 反應曲線最適法(Reaction Curve Fitting)可 以得到一階帶遲延轉移函數的參數:時延 (dead time, t0)、時間常數(time constant,

τ

)

及穩態增益(steady state gain, K)。應用上 述方法分別對 W&B 製程之2 2× 系統作轉 移函數識別。Wood & Berry 製程系統之 轉移函數為: 3 1 1 7 3 2 2 12.8 18.9 ( ) 16.7 1 21.0 1 ( ) ( ) 6.6 19.4 ( ) 10.9 1 14.4 1 s s s s e e y s s s u s y s e e u s s s − − − −  −      + +   =          + +    (12) 由於每次僅能計算單一輸出所以分別對 兩輸出作兩次2 1× 識別運算。首先對 y1 及兩個輸入作系統識別並設定三個參數 如下:(1)取樣 1000 個樣本(N=1000,取樣 間隔 1min),(2)脈衝應答係數長度 300 (n=300),(3)PLSR 潛變數(Latent Variables, lv)個數 40(lv=40)。所得結果整理如下所 示:,得到其預測殘差圖如圖二。並可以 計算所建立的模型中 40 個潛變數(lv)分別 解釋之 X 及 Y 變異及其累計值分別如表 一。並以新的常態數據(N(0,0.5), seed(3,5)) 加以驗證,如圖三。 圖二、W&B 製程中 y1建模數據之預測殘 差圖 圖三、W&B 製程中 y1建模數據之預測驗 證圖 且針對此例脈衝回歸係數(以

θ

表示) 和 PLSR1 回歸係數的關係如下:

(6)

[

b1a(1: 300) b1a(301: 600)

]

300 2

θ = × (13)

可由所得之脈衝應答係數累加計算得到 輸出 y1分別對輸入 u1及 u2的一階帶遲延

(First Order Plus Dead Time, FOPDT)模型 之一階應答如圖四及圖五。 圖四、W&B 製程 y1對 u1之階梯應答 圖五、W&B 製程 y1對 u2之階梯應答 同理,可以得到對輸出 y2之階梯應答 (參數設定和 y1相同)。將上述階梯應答曲 線所得結果可由反應曲線最適法可以得 到一階帶遲延轉移函數(14)式。 0.5225 2.5 1 6.55 2.51 2 12.79 18.91 ˆ ( ) 16.67 1 21.01 1 ˆ ( ) 6.59 19.42 10.82 1 14.42 1 s s s s e e Gp s s s e e Gp s s s − − − −       + +   =    −      + +   (14)

四、

小波分析於錯誤診斷的應用

1. 小波分析(Wavelet Analysis) 若已知母波(Mother Wavelet)及其他 所有的基底(basis)性質,則任何訊號一為 訊 號 x t( ) 透 過 尺 度 函 數 (Scaling Function)

φ

的平移(shifted versions)及小波 函數(Wavelet

表一、測試一之潛變數(lv)分別解釋之 X 及 Y 變異及其累計值

Percent Variance Captured by PLS Model

LV# 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000 21.0000 22.0000 23.0000 24.0000 25.0000 26.0000 27.0000 28.0000 29.0000 30.0000 31.0000 32.0000 33.0000 34.0000 35.0000 36.0000 37.0000 38.0000 39.0000 40.0000 ----X-Block--- This LV Total 0.2677 0.2677 0.3603 0.6281 0.4207 1.0488 0.5271 1.5759 0.4551 2.0310 0.4350 2.4660 0.4039 2.8700 0.4373 3.3072 0.3749 3.6821 0.3976 4.0797 0.4602 4.5399 0.5251 5.0651 0.4000 5.4650 0.2962 5.7613 0.3158 6.0770 0.3139 6.3909 0.3812 6.7721 0.4365 7.2087 0.4076 7.6162 0.3753 7.9915 0.2929 8.2844 0.3350 8.6193 0.2971 8.9165 0.2886 9.2051 0.3058 9.5108 0.2784 9.7892 0.4048 10.1940 0.3196 10.5136 0.2729 10.7865 0.3007 11.0872 0.2802 11.3673 0.3334 11.7008 0.3150 12.0157 0.3553 12.3711 0.2792 12.6502 0.3007 12.9510 0.3062 13.2572 0.3009 13.5581 0.2917 13.8498 0.3069 14.1567 ----Y-Block--- This LV Total 74.4144 74.4144 16.2150 90.6294 5.6221 96.2514 1.4844 97.7359 0.7935 98.5294 0.4777 99.0071 0.3042 99.3112 0.1922 99.5034 0.2044 99.7078 0.1003 99.8081 0.0479 99.8561 0.0406 99.8966 0.0231 99.9197 0.0255 99.9452 0.0161 99.9613 0.0125 99.9738 0.0078 99.9816 0.0044 99.9860 0.0033 99.9893 0.0020 99.9913 0.0016 99.9929 0.0011 99.9940 0.0008 99.9948 0.0009 99.9957 0.0008 99.9965 0.0006 99.9971 0.0006 99.9977 0.0007 99.9984 0.0005 99.9989 0.0003 99.9992 0.0002 99.9994 0.0001 99.9995 0.0001 99.9996 0.0001 99.9997 0.0001 99.9998 0.0000 99.9998 0.0000 99.9998 0.0000 99.9999 0.0000 99.9999 0.0000 99.9999

Function)ψ的解析(shifted and dilated version),可以分解並重組表示如下: 0 0 0 , , ( ) ( ) ( ) J J K J K J K J K K J K x t u φ t X ψ t =−∞ =

+

∑ ∑

(15) 其中小波函數(Wavelet Function)

ψ

及尺度 函數(Scaling Function)

φ

分別依特殊的設 計可表示如下: / 2 , ( ) 2 (2 ) J J J K t t K ψ = − ψ − (16) / 2 , ( ) 2 (2 ) J J J K t t K φ = − φ − (17) 其中J K, ∈Z,稱之為元素(atoms)。如上 構成小波解析中正交的基底(Orthonormal basis),則所有L R2( )空間中函數皆可依上 述正交元素的線性組合展開表示。 0 * , ( ) , ( ) J K J K u =

x t ⋅φ t dt (18)

(7)

* , ( ) , ( ) J K J K X =

x t ⋅Ψ t dt (19)

五、

預測式錯誤診斷與小波分析

首先,預測式錯誤診斷必須先建立製 程模型,使能執行預測的功能。本文在第 四章論及以 DPLS1 作動態系統的因果模 型(轉移函數)識別,且[測試二]已經針對本 文所討論的模擬 W&B 製程完成模型識別 的部分。因此,以下討論皆假設製程模行 為已知,且配合離散式小波轉換技術及選 取的特定母波函數(db5 及 dmey)小波執 行訊號的分解。最後,利用 PLS2 方法進 行錯誤診斷的研究。 輸出預測值 ˆy可以透過製程識別所 建立的轉移函數計算如下: 11 12 21 22 ˆ ( ) ( ) ˆ( ) ˆ ( ) ( ) D B x s G G R s y s x s G G S s       =  =         (25) 且透過小波轉換,所有輸入訊號可被分解 為低頻成分(low frequency component, l)、 中頻(middle frequency component, m)、高 頻(high frequency component, h)。可表示如 下: , l m h l m h R= +R R +R S= +S S +S (26) 由於小波重建的可加性,因此輸出預測可 由小波轉換後所得不同頻帶之訊號,分別 通過轉移函數模型之成分計算而得。其公 式計算如下: 1 11 1 2 11 2 3 11 3 4 21 1 5 21 2 6 21 3 1 21 1 2 21 2 3 21 3 4 22 1 5 22 2 6 22 3 ˆ ˆ D B X a G R a G R a G R a G S a G S a G S X b G R b G R b G R b G S b G S b G S = ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ = ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ (27) 其中 Ri、Si其中之下標 i=1~3,分別代表 由低頻至高頻訊號。常態下,理想情況的 輸出實際值與預測值相同,可表示如下: ˆ ˆ D D B B X X X X  =   =  (28) 因此,吾人可以列式如下式並執行偏 回歸分析如下: 1 11 1 2 11 2 3 11 3 4 21 1 5 21 2 6 21 3 1 21 1 2 21 2 3 21 3 4 22 1 5 22 2 6 22 3 D B X a G R a G R a G R a G S a G S a G S X b G R b G R b G R b G S b G S b G S = ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ = ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ + ⋅ ⋅ (29) 首先,在常態下建立預測變數 T2及 Q(或 SPE)正常超作上下限及輸出應答之 SPE 上下限,並分別分解各個頻帶輸入值成分 之大小及計算其正常上下限值。希望結合 小 波 分 析 的 工 具 及 貢 獻 圖 (Contribution plot)的觀念,診斷下列錯誤的發生: (1) 輸出值正常下 輸入值異常之來源量測變數及其頻 帶(低頻偏移或高頻雜訊)。 (2) 輸入值正常而輸出值異常下 輸出值異常之來源量測變數,及其訊 號的型態是否源於製程的動態異常,或源 於輸入(未能測量之)異常的擾動,而造成 輸出異常。 本文針對多變數線性系統-Wood & Berry 製程為例,先以 db5 母波( Mother Wavelet )模擬各種錯誤的預測式錯誤診斷 結果如下: 1. 各類錯誤範例測試 執行 PLS 之前,設計預測變數 X 矩 陣及被預測變數 Y 矩陣如下:

[

]

11 1 11 3 12 12 3 21 1 11 3 22 12 3 n12 X g r g r g r g r g r g r g r g r × = L L L L L (30) D B X Y X   =   (31)

(8)

其中g rij k為輸入 R 之第 k(k=1~3)頻帶成分 k r 經過gij轉移函數(i, j=1~2)的預測,所得 到 的 小 波 轉 換 輸入 (Wavelet transformed input data)訊號對應之輸出預測值,包含 12 個成分(變數),n 為取樣個數。而被預 測變數 Y 為輸出實際量測值,透過回歸之 線性分析可以了解預測變數 X 及被預測 變數 Y 是否皆處於常態,亦可用以偵測錯 誤的來源。 將 X 及 Y 經過均值移除(mean center) 與正規化的前處理後,以 NIPALS 演算法 建立 PLS 模型,得到 X 的權重(weights) w 及其負荷(loadings) p,同時得到相對應的 Y 負荷(loadings) q 及內相關係數(inner relation coefficients) b。並可了解 X 空間解 釋 Y 空間的變異量之比例大小(ssq)。經過 上述模型的建立,任何新的樣本都可以上 列參數轉換至 X 變量 t 及 Y 變量 u。 對於前變數個數的選擇,實際製程中 由於大量變數產生的線性相關問題,必須 將變異極小的變量移除。變異極小的變量 除了可能代表無意義的雜訊外,計算上可 以避免 T2 統計量運算時因分母過小造成 計算上的困難,例如導致 T2統計量於正常 操作下反而異常變大的情況(T. Kourti and J. F. MacGregor, 1995)[2]。但是 P. Geladi and B. Kowalski (1986)[3]曾經提及移除變 量使訊息失真的可能性。本文沒有刪除任 何變量(即變數個數 lv=12)。並將結果整理 用以偵錯。 (1) 正常操作狀態(樣本數 n=400) (輸入值 R, S~N(0,0.5),其輸入雜訊 (input noise)分別皆為 N(0,0.01),而輸出雜 訊為(Output noise) N(0,0.2))。且其對常態 0 50 100 150 200 250 300 350 400 0 10 20 30 T

2 and Q chart for X-block and Y-block T 2 of X-b lo c k Tapha 0 50 100 150 200 250 300 350 400 0 1 2 3x 10 -29 Q o f X-b lo c k 10,7,3*std(Qn) Qapha 0 50 100 150 200 250 300 350 400 0 0.1 0.2 0.3 0.4 Q of o ut p ut sample number 10,7,3*std(Qn) Qyapha 圖六、常態數據 T2及 Q 管制圖 0 100 200 300 400 -3 -2 -1 0 1 2 3

Y1: Actual value(-) vs.Predicted Value (--)

Ou t p u t 0 100 200 300 400 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 sample number Y1 Prediction Error Re sid u al 0 100 200 300 400 -3 -2 -1 0 1 2 3

Y2 :Actual value(-) vs.Predicted Value (--)

0 100 200 300 400 -0.6 -0.4 -0.2 0 0.2 0.4 Y2 Prediction Error sample number 圖七、常態數據 Y 殘差管制圖 0 100 200 300 400 -2 0 2

Contribution plot for output Y1

g r11 1 0 100 200 300 400 -2 0 2 g r11 2 0 100 200 300 400 -1 0 1 g r11 3 0 100 200 300 400 -2 0 2 g s21 1 0 100 200 300 400 -2 0 2 g s21 2 0 100 200 300 400 -2 0 2 g s21 3 0 100 200 300 400 -0.2 0 0.2 g21 r1 0 100 200 300 400 -0.1 0 0.1 g21 r2 0 100 200 300 400 -0.05 0 0.05 g r21 3 0 100 200 300 400 -0.5 0 0.5 g22 s 1 0 100 200 300 400 -0.5 0 0.5 g22 s 2 0 100 200 300 400 -1 0 1 g s22 3 圖八、常態數據各頻帶對 Y1 貢獻管制圖

(9)

0 100 200 300 400 -0.2

0 0.2

Contribution plot for output Y2

g r11 1 0 100 200 300 400 -0.05 0 0.05 g r11 2 0 100 200 300 400 -0.05 0 0.05 g r11 3 0 100 200 300 400 -0.5 0 0.5 g s21 1 0 100 200 300 400 -1 0 1 g s21 2 0 100 200 300 400 -2 0 2 g s21 3 0 100 200 300 400 -1 0 1 g r21 1 0 100 200 300 400 -1 0 1 g r21 2 0 100 200 300 400 -0.5 0 0.5 g r21 3 0 100 200 300 400 -1 0 1 g s22 1 0 100 200 300 400 -2 0 2 g s22 2 0 100 200 300 400 -0.5 0 0.5 g s22 3 圖九、常態數據各頻帶對 Y2 貢獻管制圖 樣本作變量轉換後,所得 T2 及 Q 管制圖 如圖六。(以下設定 T2管制圖的上限值取 99%信賴區間,alpha=0.01)。為了簡約篇 幅,以下僅示範一種因輸入量測偏差引起 的異常狀況。 給定第一個輸入值 R 發生於 t=100(時 間獲取樣點)的量測偏移錯誤(5 個單位), 意即輸入之偏移錯誤未通過製程,僅影響 預測值。 0 50 100 150 200 250 300 350 400 0 20 40 60 T

2 and Q chart for X-block and Y-block

T 2 o f X-b lo c k Ta 0 50 100 150 200 250 300 350 400 0 1 2x 10 -28 Q o f X-b lo c k 10,7,3*std(Qn),Qa 0 50 100 150 200 250 300 350 400 0 10 20 30 40 Q o f o ut p ut sample number 10,7,3*std(Qn), Qya 圖十、輸入量測偏差異常數據 T2及 Q 管 制圖 0 100 200 300 400 -4 -2 0 2 4 6

8Y1: Actual value(-) vs.Predicted Value (--)

Ou t p ut 0 100 200 300 400 -5 -4 -3 -2 -1 0 1 sample number Y1 Prediction Error R e sid u al 0 100 200 300 400 -4 -2 0 2 4

6Y2 :Actual value(-) vs.Predicted Value (--)

0 100 200 300 400 -3 -2 -1 0 1 Y2 Prediction Error sample number 圖十一、輸入偵測器偏差異常數據 Y 殘差 管制圖 0 100 200 300 400 -5 0 5

Contribution plot for output Y1

g r111 0 100 200 300 400 -2 0 2 g r112 0 100 200 300 400 -1 0 1 g1 1 r3 0 100 200 300 400 -2 0 2 g s21 1 0 100 200 300 400 -2 0 2 g s21 2 0 100 200 300 400 -2 0 2 g s21 3 0 100 200 300 400 -0.5 0 0.5 g r211 0 100 200 300 400 -0.1 0 0.1 g r212 0 100 200 300 400 -0.05 0 0.05 g2 1 r3 0 100 200 300 400 -0.5 0 0.5 g s22 1 0 100 200 300 400 -0.5 0 0.5 g s22 2 0 100 200 300 400 -1 0 1 g s22 3 圖十二、輸入偵測器偏差異常各頻帶對 Y1 貢獻管制圖 由第一個管制圖(圖十)顯示輸入有 發生錯誤,且輸出的 SPE 圖超出上限。 0 100 200 300 400 -1 0 1

Contribution plot for output Y2

g r111 0 100 200 300 400 -0.05 0 0.05 g r112 0 100 200 300 400 -0.1 0 0.1 g11 r3 0 100 200 300 400 -0.5 0 0.5 g21 s1 0 100 200 300 400 -1 0 1 g s21 2 0 100 200 300 400 -2 0 2 g s21 3 0 100 200 300 400 -5 0 5 g r211 0 100 200 300 400 -1 0 1 g r212 0 100 200 300 400 -0.5 0 0.5 g21 r3 0 100 200 300 400 -1 0 1 g s22 1 0 100 200 300 400 -2 0 2 g s22 2 0 100 200 300 400 -0.5 0 0.5 g22 s3 圖十三、輸入偵測器偏差異常各頻帶對 Y2 貢獻管制圖 而第二個管制圖(圖十一)顯示輸出之實 際值皆在上下限的範圍內,但預測值卻受 到輸入值量測錯誤的影響而使輸出預測 值分布在正常範圍之外。然而,輸入錯誤 的來源可由第三及第四個管制圖(圖十二

(10)

及十三)作診斷。由於此錯誤為低頻的偏 移,因此最後兩個管制圖中和 r1有關之變 數皆發生異常偏移正常範圍。由於來自輸 入的量測錯誤會影響到所有和此量測量 相關的所有輸出之預測,所以錯誤會成對 出現。此乃本類錯誤的重要參考特徵。 (3)輸入的高頻量測錯誤(N(0,10)) 0 50 100 150 200 250 300 350 400 0 20 40 60 80 T

2 and Q chart for X-block and Y-block T 2 of X-b lo c k Tapha 0 50 100 150 200 250 300 350 400 0 1 2x 10 -28 Q o f X-b lo c k 10,7,3*std(Qn) Qapha 0 50 100 150 200 250 300 350 400 0 5 10 Q o f o u t p u t sample number 10,7,3*std(Qn) Qyapha 圖十四、輸入偵測器高頻異常數據 T2及 Q 管制圖 給定第一個輸入值 R,一發生於 t=100(時間獲取樣點)的高頻量測錯誤(變 異為 10 單位(N(0,10))的高頻成分),意即 輸入之偏移錯誤未通過製程,僅影響預測 值。 0 100 200 300 400 -4 -2 0 2

4Y1: Actual value(-) vs.Predicted Value (--)

O ut p ut 0 100 200 300 400 -3 -2 -1 0 1 2 3 sample number Y1 Prediction Error R esid u al 0 100 200 300 400 -3 -2 -1 0 1 2

3Y2 :Actual value(-) vs.Predicted Value (--)

0 100 200 300 400 -2 -1 0 1 2 Y2 Prediction Error sample number 圖十五、輸入偵測器高頻異常數據 Y 殘差 管制圖 0 100 200 300 400 -2 0 2

Contribution plot for output Y1

g r111 0 100 200 300 400 -2 0 2 g r112 0 100 200 300 400 -5 0 5 g r11 3 0 100 200 300 400 -2 0 2 g s21 1 0 100 200 300 400 -2 0 2 g s21 2 0 100 200 300 400 -2 0 2 g s21 3 0 100 200 300 400 -0.2 0 0.2 g r211 0 100 200 300 400 -0.1 0 0.1 g r212 0 100 200 300 400 -0.5 0 0.5 g r21 3 0 100 200 300 400 -0.5 0 0.5 g s22 1 0 100 200 300 400 -0.5 0 0.5 g s22 2 0 100 200 300 400 -1 0 1 g s22 3 圖十六、輸入偵測器高頻異常對 Y1 貢獻 管制圖 0 100 200 300 400 -0.2 0 0.2

Contribution plot for output Y2

g11 r1 0 100 200 300 400 -0.05 0 0.05 g r112 0 100 200 300 400 -0.5 0 0.5 g11 r3 0 100 200 300 400 -0.5 0 0.5 g2 1 s 1 0 100 200 300 400 -1 0 1 g s21 2 0 100 200 300 400 -2 0 2 g s21 3 0 100 200 300 400 -1 0 1 g r211 0 100 200 300 400 -1 0 1 g r212 0 100 200 300 400 -5 0 5 g r21 3 0 100 200 300 400 -1 0 1 g2 2 s 1 0 100 200 300 400 -2 0 2 g s22 2 0 100 200 300 400 -0.5 0 0.5 g s22 3 圖十七、輸入偵測器高頻異常對 Y2 貢獻 管制圖 由第一個管制圖(圖十四)顯示輸入 有發生錯誤,且輸出的 SPE 圖超出上限。 而第二個管制圖顯示輸出之實際值皆在 上下限的範圍內,但預測值卻受到輸入值 量測錯誤的影響而使輸出預測值高頻的 震盪且分布在正常範圍之外。然而,輸入 錯誤的來源可由第三及第四個管制圖作 診斷。由於此錯誤為高頻雜訊,因此最後 兩個管制圖中和 r3有關之變數皆發生異常 偏移正常範圍。由於來自輸入的量測錯誤 會影響到所有和此量測量相關的所有輸 出之預測,所以錯誤會成對出現。此乃本 類錯誤的重要參考特徵。 [比較] 輸入的高頻量測錯誤(N(0,1)) 給 定 第 一 個 輸 入 值 R , 一 發 生 於 t=100(時間或取樣點)的高頻量測錯誤(變 異為 1 單位(N(0,1))的高頻成分),意即輸

(11)

入之偏移錯誤未通過製程,僅影響預測 值。 0 50 100 150 200 250 300 350 400 0 20 40

60 T2 and Q chart for X block and Y block of X-T h 0 50 100 150 200 250 300 350 400 0 1 2 3 4x 10 Q of 10 7 3*std(Q ) Qh 0 50 100 150 200 250 300 350 400 0 0 5 1 Q of sample number 10 7 3*std(Q ) Q h 圖十八、輸入的高頻量測錯誤(N(0,1))T2 及 Q 管制圖 0 100 200 300 400 -3 -2 -1 0 1 2

3Y1: Actual value(-) vs.Predicted Value (--)

O ut p ut 0 100 200 300 400 -1 -0.5 0 0.5 1 sample number Y1 Prediction Error R e sid u al 0 100 200 300 400 -3 -2 -1 0 1 2

3Y2 :Actual value(-) vs.Predicted Value (--)

0 100 200 300 400 -1 -0.5 0 0.5 1 Y2 Prediction Error sample number 圖十九、輸入的高頻量測錯誤(N(0,1))Y 殘差管制圖 由於此雜訊變異較上述小 10 倍,所 以第一個管制圖(圖十八)輸入 X 的 T2 Q 超出上限的分界並不明顯,且輸出的 SPE 圖超出上限雖可判斷,但起始點亦不 0 100 200 300 400 -2 0 2

Contribution plot for output Y1

g r111 0 100 200 300 400 -2 0 2 g r11 2 0 100 200 300 400 -1 0 1 g r113 0 100 200 300 400 -2 0 2 g s21 1 0 100 200 300 400 -2 0 2 g s21 2 0 100 200 300 400 -2 0 2 g s21 3 0 100 200 300 400 -0.2 0 0.2 g r211 0 100 200 300 400 -0.1 0 0.1 g r21 2 0 100 200 300 400 -0.1 0 0.1 g r213 0 100 200 300 400 -0.5 0 0.5 g s22 1 0 100 200 300 400 -0.5 0 0.5 g s22 2 0 100 200 300 400 -1 0 1 g s22 3 圖二十、輸入的高頻量測錯誤(N(0,1))輸 入對輸出 Y1 貢獻管制圖 0 100 200 300 400 -0.2 0 0.2

Contribution plot for output Y2 g r111 0 100 200 300 400 -0.05 0 0.05 g11 r2 0 100 200 300 400 -0.1 0 0.1 g r11 3 0 100 200 300 400 -0.5 0 0.5 g s21 1 0 100 200 300 400 -1 0 1 g s21 2 0 100 200 300 400 -2 0 2 g s213 0 100 200 300 400 -1 0 1 g21 r1 0 100 200 300 400 -1 0 1 g21 r2 0 100 200 300 400 -1 0 1 g21 r3 0 100 200 300 400 -1 0 1 g22 s 1 0 100 200 300 400 -2 0 2 g22 s 2 0 100 200 300 400 -0.5 0 0.5 g s223 圖二十一、輸入的高頻量測錯誤(N(0,1)) 輸入對輸出 Y2 貢獻管制圖 顯著。而第二個管制圖顯示輸出之實際值 與預測值皆在上下限的範圍內,但兩量測 輸出各自的預測殘差值卻高頻地超出正 常範圍,但不及上例明顯。然而,輸入錯 誤的來源可由第三及第四個管制圖作診 斷。由於此錯誤為高頻雜訊,因此最後兩 個管制圖中和 r3有關之變數皆應發生異常 偏移正常範圍,但十分不顯著。可見其偵 查的敏銳度(sensitive)亦有其限制存在。

六、 結論

本文結合傳統 FDI 方法的概念,用 較為簡化的方法-動態 PLSR1 (dynamic PLSR 1),建 構製程的因 果模型 (Causal model);配合多變量統計原理為基礎的知 識 庫 方 法 (Knowledge-based method)- PLSR2,執行偵錯的功能。輔以離散式小 波轉換(Discrete Wavelet Transform, DWT) 技術,依頻率變化將訊號解析為趨勢變化 (trend component) 、 週 期 變 化 (seasonal component)及高頻雜訊(stationary random noise component)三種分類成分。綜合上 述,將製程輸入(Input)數據透過小波分析 及模式預測,得到輸入訊號各頻帶對輸出 的貢獻及輸出預測值(Predicted output) ; 並與製程輸出值作比較,作為錯誤診斷的 基準及法則(Fault diagnosis rules)。因此, 稱之為預測式錯誤診斷(Predictive Fault Diagnosis) 。 最 後 , 以 多 變 數 線 性 系 統 -Wood & Berry 製程為例,測試其於各種 異常狀況下的診斷效果。結果顯示,可直

(12)

接由製程輸入資料預測相對輸出錯誤發 生的時間,並探討訊號趨勢偏移(bias or drift)及高頻擾動的發生,以診斷並推論造 成錯誤的原因。

此種方法有別於前述的知識庫方

,不需要針對所有已知且可能發生的錯 誤狀態收集樣本並建模及辨識

。本法

直接 由常態建模及簡單的法則,便可作正確的 推論診斷

。但是實際應用於工廠的錯誤

偵 錯 及 診 斷 而 言

, 因 果 模 型(Causal model)

數據的取得十分困難。然而

, PLSR1

製程識別的步驟解決了系統動

態的問題

,輸入各頻帶訊號透過轉移函數 的轉換,而與系統輸出完全線性相關

。至

於反應器內部錯誤的診斷

,是後續值得 探討的問題

七、參

[1] T. Kourti and J. F. MacGregor, Process analysis, monitoring and diagnosis, using multivariate projection methods, Chemometrics and Intelligent Laboratory Systems, 28 (1995) 3-21.

[2] W. Ku et al., Disturbance detection and isolation by dynamic principal component analysis, Chemometrics and Intelligent Laboratory Systems, 30 (1995) 179-196. [3] P. Geladi and B. Kowalski, Partial Least-Squares Regression: A tutorial, Analytica Chimica Acta, 185 (1986) 1-17. [4] S. Wold et al., PLS-regression: a basic tool of chemometrics, Chemometrics and Intelligent Laboratory Systems, 58 (2001) 109-130.

[5] H. Wold, Estimation of Principle Components and Related Models by Iterative Least Squares, Multivariate Analysis, R. P. Krishnaiah, ed., Academic Press. New York (1966), p.391.

[6] R. Manne, Analysis of Two Partial Least Squares Algorithms for Multivariate Calibration, Chemometrics and Intelligent Laboratory Systems 2(1987) 283.

[7] Inge. S. Helland, Some theoretical aspects of partial least squares regression, Chemometrics and Intelligent Laboratory

Systems 58(2001) 97-107.

[8] N. L. Ricker, The use of biased least-squares estimators for parameters in discrete-time pulse response models, Ind. Engng Chem. Res. 27, 343-350.

[9] M. H. Kaspar, Model identification for chemical process control, Ph.D. thesis, University of Wisconsin, Madison (1992). [10] M. H. Kaspar and W. H. Ray ,Chemometric methods for process monitoring and high performance controller design, A.I.Ch.E. J. 38 (1992) 1593-1608. [11] S. J. Qin and T. J. McAvoy, , Nonlinear PLS Modeling using Neural Networks, Comput. Chem. Eng. 16 (4) 379 (1992b). [12] U. Kruger et. Al., An Alternative PLS Algorithm for the Monitoring of Industrial Process, Proceedings of the American Control Conference, Arlington, VA June 25-27, 2001.

[13] S. J. Qin, and T. J. McAvoy, Nonlinear PLS Modeling Via A Neural Net PLS Approach, Comput. Chem. Eng. 20 (2) 147-159 (1996).

[14] L. Ljung, Issues in system identification, IEEE Control Sysytem Magazine, 11 (1), 25-29 (1991).

[15] J. F. Macgregor and T. Kourti, Statistical process control of multivariate processes, Control Eng. Practice, vol. 3, No. 3, pp. 403-414 (1995).

[16] A. Bakhtazad et al., Process data de-noising using wavelet transform, Intell. Data Analysis 3 (1999) 267-285.

[17] K. A. Kosanovich and M. J. Piovoso, PCA of Wavelet Transformed Process Data for Monitoring, Intell. Data Analysis 1 (1997) 85-99.

[18] B. M. Wise et al., A theorical basis for the use of principlal component models for monitoring Multivariate process, Process Control and Quality, (1990) 1, pp.41-51.

[19] 吳敏綺,應用次空間程序識別做動態

系統之監控及偵錯,碩士論文,國立台灣 大學化學工程學研究所,台北 (2002)。

參考文獻

相關文件

Textbook Chapter 33.4 – Finding the closest pair of points.. Closest Pair of

‹ ‹ A product term A product term implicant implicant is called a prime is called a prime implicant implicant if it cannot be combined with another term to. if it cannot be

  本作品希望透過 數位科技 與 藝術創作 產生的激盪,讓人、環境、科技 之間形成更和諧、優質的生活形態。..   因此將利用數位科技的 互動性

Reading: Stankovic, et al., “Implications of Classical Scheduling Results for Real-Time Systems,” IEEE Computer, June 1995, pp.. Copyright: All rights reserved, Prof. Stankovic,

55 統計數據在數據庫作出修訂後,並沒同步修訂文字版的數據內容,對使 用者造成混淆錯引數字(例如數據庫中顯示 2015 年會議及展覽數目有 1263 項,惟文字版仍寫著 909

He proposed a fixed point algorithm and a gradient projection method with constant step size based on the dual formulation of total variation.. These two algorithms soon became

背景:一名小學生家長投訴學校在沒有通 知家長的情況下,向網絡程式供應商提供

mathematical statistics, statistical methods, regression, survival data analysis, categorical data analysis, multivariate statistical methods, experimental design.