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