KPCAコード

KPCAアルゴリズムのTE過程での故障診断の応用

KPCAアルゴリズム紹介

KPCAコード

データ前処理

核行列を計算する

.

中心化コア行列

共分散行列の固有値分解

により、主元の個数

が決定される.

は、特徴ベクトルを特徴値の大きさ順に並べ替える

.

単位化特徴ベクトル

試験データを再構築する

制御限の設定

故障データの検出

描画

.検出率を求める

.
KPCAアルゴリズムの紹介

    KPCA          ，                。      ：                                              ，    PCA                ，                       。
             ：
    1.      
    2.    
    3.            
    4.          
    5.       
    6.       SPE    T2   
    7.          
    8.      
    9.      
    10.             
    11.       
    12.       SPE    T2   
    
                 KPCA   TE         ，    ：

KPCAコード
データ前処理

patterns=zscore(train);            
%%            
m=mean(train);
s=std(train);
n=size(test,1);%            
test_patterns=(test-repmat(m,n,1))./repmat(s,n,1);
train_num=size(patterns,1);%    （    ）
test_num=size(test_patterns,1);%    （    ）
cov_size = train_num;
cov_size2 = test_num;

核行列の計算

for i=1:cov_size
for j=i:cov_size
K(i,j) = exp(-norm(patterns(i,:)-patterns(j,:))^2/rbf_var);
K(j,i) = K(i,j);
end
end
unit = ones(cov_size, cov_size)/cov_size;

中心化核行列

K_n = K - unit*K - K*unit + unit*K*unit;

共分散行列の固有値分解

[evectors,evaltures] = eig(K_n);
[x,index]=sort(real(diag(evaltures))); 
 evals=flipud(x);%    
 index=flipud(index);
evaltures=real(diag(evaltures));

主要素の個数を確定します

umpc=1;
 while sum(evals(1:numpc))/sum(evals)<0.85
    numpc=numpc+1;
 end
 a=numpc;

特徴ベクトルを特徴値の大きさ順に並べ替えます.

evectors=evectors(:,index);

単位化特徴ベクトル

for i=1:cov_size,
evecs(:,i) = evectors(:,i)/sqrt(evals(i));
end 
index=numpc;
train_kpca=K_n * evecs(:,1:index(1)); %%          
scoresp=K_n*evecs; %%%                      

試験データの再構築

unit_test = ones(test_num,cov_size)/cov_size;
K_test = zeros(test_num,cov_size);
for i=1:test_num
    for j=1:cov_size
        K_test(i,j) = exp(-norm(test_patterns(i,:)-patterns(j,:))^2/rbf_var);
    end
end
K_test_n = K_test - unit_test*K - K_test*unit + unit_test*K*unit;%        
test_kpca = K_test_n * evecs(:,1:index(1));%%         
scoresp2=K_test_n*evecs; %%%                      

制御限の設定

for i=1:cov_size
        T12(i)=(K_n(i,:)*evecs(:,1:a)*inv(diag(evals(1:index(1)))/train_num)*evecs(:,1:a)'*K_n(i,:)');%      
end
Q1=(sum(scoresp.^2, 2)-sum(train_kpca.^2, 2))';
T2cfdLimit = ksdensity(T12,0.99,'function','icdf');%icdf：inverse cdf         
SPEcfdLimit = ksdensity(Q1,0.99,'function','icdf');

故障データの検出

for i=1:cov_size2
        T22(i)=(K_test_n(i,:)*evecs(:,1:a)*inv(diag(evals(1:index(1)))/train_num)*evecs(:,1:a)'*K_test_n(i,:)');%      
        Q2=(sum(scoresp2.^2, 2)-sum(test_kpca.^2, 2))';

図形描画

figure;
subplot(2,1,1);
plot(1:cov_size,T12,'blue');
xlabel('     (  )');
ylabel('T2');
hold on;
line([0,cov_size],[T2cfdLimit,T2cfdLimit],'LineStyle','- -','color','red');
legend('   ','  ');

subplot(2,1,2);
plot(1:cov_size,Q1,'blue');
xlabel('     （  ）');
ylabel('Q');
hold on;
line([0,cov_size],[SPEcfdLimit,SPEcfdLimit],'lineStyle','- -','Color','red');
legend('   ','  ');

figure;
subplot(2,1,1);
plot(1:cov_size2,T22,'blue');
xlabel('   (  )');
ylabel('T2');
hold on;
line([0,cov_size2],[T2cfdLimit,T2cfdLimit],'LineStyle','- -','color','red');
legend('   ','  ');

検出率を求める

cout_T=0;
cout_SPE=0;
cout_T1=0;
cout_SPE1=0;

for i=1:cov_size2
  if i<=160 && T22(i)>T2cfdLimit
     cout_T1=cout_T1+1;
  end 
  if i<=160 && Q2(i)>SPEcfdLimit
     cout_SPE1=cout_SPE1+1;
  end
  if i>160 && T22(i)>T2cfdLimit
     cout_T=cout_T+1;
  end
  if i>160 && Q2(i)>SPEcfdLimit
     cout_SPE=cout_SPE+1;
  end
end
cout_T_F=cout_T1/160*100
cout_T=cout_T/(cov_size2-160)*100
cout_SPE_F=cout_SPE1/160*100
cout_SPE=cout_SPE/(cov_size2-160)*100