PCA故障診斷步驟

1、基于PCA算法的故障診斷步驟離線PCA監(jiān)測(cè)模型的計(jì)算步驟：（1） 選擇監(jiān)控變量，收集正常工況下的各變量的樣本，記為訓(xùn)練樣本數(shù)據(jù)X_train和檢驗(yàn)數(shù)據(jù)X_test；X_train =x11 x12x1mx21 x22x2mxn1 xn2xnmX_train為n×m矩陣，即n個(gè)樣本，m個(gè)觀測(cè)變量（即以列向量來(lái)看的話，為一個(gè)觀測(cè)變量各個(gè)采樣點(diǎn)的值）對(duì)樣本數(shù)據(jù)X_train和檢驗(yàn)數(shù)據(jù)X_test進(jìn)行中心化和標(biāo)準(zhǔn)化處理得到和Xstest；中心化處理：按列對(duì)X_train減去觀測(cè)變量的均值xj=xj-xj 觀測(cè)變量某一采樣點(diǎn)的值減去這一觀測(cè)變量所有采樣點(diǎn)的平均值xj=1ni=1nxij 求取一

2、列（即某一觀測(cè)變量）的平均值標(biāo)準(zhǔn)化處理：對(duì)X_train除以觀測(cè)變量的標(biāo)準(zhǔn)差（按列（觀測(cè)變量）進(jìn)行）xs j=xjsjsj=i=1n（xij-xj）2n-1 標(biāo)準(zhǔn)差求出標(biāo)準(zhǔn)化矩陣的協(xié)方差矩陣；的協(xié)方差矩陣對(duì)為：=1n-1XstainT*Xstrain（2） 對(duì)進(jìn)行特征分解,求得特征值1,2,m（1>2>>m）及其對(duì)應(yīng)的特征向量p1,p2,pm（負(fù)荷向量）；（3） 確定主元個(gè)數(shù), 確定了主元個(gè)數(shù)k,就得到了k個(gè)特征值1>2>>k,及其對(duì)應(yīng)的特征向量p1,p2,pk；A:累計(jì)貢獻(xiàn)率法：前k個(gè)主元的累積方差貢獻(xiàn)率為： i=1kii=1mi當(dāng)前k個(gè)主元的累積方差貢獻(xiàn)

3、率達(dá)到85%，則主元個(gè)數(shù)取k值B:交叉檢驗(yàn)估計(jì)法：將采集到的數(shù)據(jù)分成k個(gè)部分，1部分?jǐn)?shù)據(jù)用來(lái)建立主元模型,剩下的k-1部分用來(lái)作為檢驗(yàn)數(shù)據(jù)去檢驗(yàn)所建的模型。如此,建立若干個(gè)不同主元個(gè)數(shù)的模型,并測(cè)試所建立的模型,從中選取一個(gè)通過檢驗(yàn)后誤差最小的模型的主元個(gè)數(shù)作為系統(tǒng)主元個(gè)數(shù)。（4） 建立PCA主元模型，并進(jìn)行交叉驗(yàn)證以確定誤差最小按照ti=XstrainPi,求出第i個(gè)主元，并依據(jù)Xstrain=t1p1T+t2p2T+tkpkT =Xstrainp1p1T+Xstrainp2p2T+XstrainpkpkT求出其主元模型用Xstest帶入得到另一主元模型Xstest，依據(jù)E=X-X,求出模型

4、誤差，確定模型誤差最小的那個(gè)模型即為主元模型。（5） 計(jì)算T2統(tǒng)計(jì)量控制限和SPE統(tǒng)計(jì)量控制限；對(duì)于樣本個(gè)數(shù)為n,主元個(gè)數(shù)為k的過程變量X_train, T2統(tǒng)計(jì)量服從自由度為k和n一k的F分布,則置信度為的T2統(tǒng)計(jì)量控制上限為:T2UCL=k(n-1)n-kF(k,n-k)或T2UCL=k(n-1)(n+1)n(n-k)F(k,n-k) 檢驗(yàn)水平為的SPE統(tǒng)計(jì)量控制上限為：QUCL=11+h0C221+2h0(h0-1)121h01=i=k+1mi, 2=i=k+1mi2,3=i=k+1mi3 , h0=1-223322, C是與（1-）分位點(diǎn)對(duì)應(yīng)的標(biāo)準(zhǔn)差在線過程監(jiān)測(cè)與故障診斷步驟：（1）

5、采集第i時(shí)刻的在線實(shí)時(shí)數(shù)據(jù)Xi（Xi為1×m矩陣），并進(jìn)行中心化和標(biāo)準(zhǔn)化處理得到Xis；（2） 按照t(i)(1×k)=Xis(1×m)pk(m×k),求出Xis的得分向量，依據(jù)Xis=tipkT(k×m),求出PCA模型估計(jì)量Xis,這里pk=p1k,p2k,pkk；（3） 計(jì)算Xis的T2統(tǒng)計(jì)量和SPE統(tǒng)計(jì)量，并畫出T2統(tǒng)計(jì)量和SPE統(tǒng)計(jì)量的控制圖；T(i)2=t(i)k-1t(i)T=Xispkk-1pkTXisTQ(i)=Xis(I-pkpkT)XisT（4） 將上述計(jì)算結(jié)果與T2統(tǒng)計(jì)量控制限和SPE統(tǒng)計(jì)量控制限比較,以檢測(cè)過程運(yùn)行有無(wú)

6、異常，當(dāng)有異常狀態(tài)發(fā)生時(shí),繪制貢獻(xiàn)圖，找出與故障相關(guān)的系統(tǒng)變量：1) 檢查每個(gè)觀測(cè)值x的標(biāo)準(zhǔn)化得分tii, 并確定造成失控狀態(tài)的r(r<a)個(gè)得分: tii<T2UCLa；2) 計(jì)算每個(gè)變量xi相對(duì)于失控得分ti的貢獻(xiàn)率是:conti,j=tiipi,jxj3) 當(dāng)conti,j是負(fù)時(shí),設(shè)它為零；4) 計(jì)算第j個(gè)過程變量xj的總貢獻(xiàn)率：CONTj=i=1r(conti,j)5) 把所有m個(gè)過程變量xj的CONTj畫在一個(gè)曲線圖上。PCA_TE仿真程序：%TE過程的傳統(tǒng)主元分析在Matlab中的仿真程序%建立模型：%載入模型數(shù)據(jù)，以故障11為例Xtrain = pkx101; Xtr

7、ain = double(Xtrain);%載入測(cè)試數(shù)據(jù)Xtest = pkx102;Xtest = double(Xtest);%標(biāo)準(zhǔn)化處理：X_mean = mean(Xtrain); %按列求Xtrain平均值 X_std = std(Xtrain); %求標(biāo)準(zhǔn)差 X_row,X_col = size(Xtrain); %求Xtrain行、列數(shù) % for i = 1:X_col%Xtrain(:,i)=(Xtrain(:,i) - X_mean(i)./X_std(i);%Xtest(:,i) = (Xtest(:,i) - X_mean(i)./X_std(i);% end Xtra

8、in=(Xtrain-repmat(X_mean,X_row,1)./repmat(X_std,X_row,1);%求協(xié)方差矩陣sigmaXtrain = cov(Xtrain);%對(duì)協(xié)方差矩陣進(jìn)行特征分解，lamda為特征值構(gòu)成的對(duì)角陣，T的列為單位特征向量，且與lamda中的特征值一一對(duì)應(yīng)：T,lamda = eig(sigmaXtrain); % disp('特征根（由小到大）');% disp(lamda);% disp('特征向量：');% disp(T); %取對(duì)角元素(結(jié)果為一列向量)，即lamda值，并上下反轉(zhuǎn)使其從大到小排列，主元個(gè)數(shù)初值為1，

9、若累計(jì)貢獻(xiàn)率小于90%則增加主元個(gè)數(shù)D = flipud(diag(lamda); num_pc = 1; while sum(D(1:num_pc)/sum(D) < 0.9 num_pc = num_pc +1;end %取與lamda相對(duì)應(yīng)的特征向量P = T(:,X_col-num_pc+1:X_col); %求置信度為99%時(shí)的T2統(tǒng)計(jì)控制限 T2UCL=num_pc*(X_row-1)*(X_row+1)*finv(0.99,num_pc,X_row - num_pc)/(X_row*(X_row - num_pc);%置信度為99%的Q統(tǒng)計(jì)控制限for i = 1:3 th

10、eta(i) = sum(D(num_pc+1:X_col).i);endh0 = 1 - 2*theta(1)*theta(3)/(3*theta(2)2);ca = norminv(0.99,0,1);QUCL = theta(1)*(h0*ca*sqrt(2*theta(2)/theta(1) + 1 + theta(2)*h0*(h0 - 1)/theta(1)2)(1/h0); %在線監(jiān)測(cè)：%標(biāo)準(zhǔn)化處理n = size(Xtest,1);Xtest=(Xtest-repmat(X_mean,n,1)./repmat(X_std,n,1);%求T2統(tǒng)計(jì)量，Q統(tǒng)計(jì)量r,y = size(

11、P*P');I = eye(r,y);T2 = zeros(n,1);Q = zeros(n,1);for i = 1:nT2(i)=Xtest(i,:)*P*inv(lamda(20-num_pc+1:20,20-num_pc+1:20)*P'*Xtest(i,:)' Q(i) = Xtest(i,:)*(I - P*P')*Xtest(i,:)' end%繪圖 figure subplot(2,1,1); plot(1:n,T2,'k'); title('主元分析統(tǒng)計(jì)量變化圖'); xlabel('采樣數(shù)

12、9;); ylabel('T2'); hold on; line(0,n,T2UCL,T2UCL,'LineStyle','-','Color','r'); subplot(2,1,2); plot(1:n,Q,'k'); xlabel('采樣數(shù)'); ylabel('Q'); hold on; line(0,n,QUCL,QUCL,'LineStyle','-','Color','r'); %貢獻(xiàn)圖%1

13、.確定造成失控狀態(tài)的得分S = Xtest(400,:)*P(:,1:num_pc);r = ;for i = 1:num_pc if S(i)2/lamda(i) > T2UCL/num_pc r = cat(2,r,i); endend%2.計(jì)算每個(gè)變量相對(duì)于上述失控得分的貢獻(xiàn)cont = zeros(length(r),20);for i = length(r) for j = 1:20 cont(i,j) = abs(S(i)/D(i)*P(j,i)*Xtest(400,j); endend%3.計(jì)算每個(gè)變量的總貢獻(xiàn)CONTJ = zeros(20,1);for j = 1:20 CONTJ(j) = sum(cont(:,j);end%4.計(jì)算每個(gè)變量對(duì)Q的貢獻(xiàn)e = Xtest(400,:)*(I - P*P')