KPCA故障诊断matlab实现

%正常模型建立 Xtrain=dlmread('d00_te.dat'); %取训练数据阵的行和列数 [XX_row,XX_col]=size(Xtrain); %平均值 Xt_mean = mean(Xtrain); %方差 Xt_std = std(Xtrain); %训练集矩阵每一列减去平均值除以方差 for i = 1:XX_col Xtrain(:,i) = (Xtrain(:,i)-Xt_mean(i))./Xt_std(i); end %计算训练数据的核矩阵K for i=1:XX_row for j=1:XX_row kX(i,j)=exp(-(norm(Xtrain(i,:)-Xtrain(j,:))).^2/400); end end lmX=1/XX_row*ones(XX_row);%定义训练数据的IN %对训练数据的核矩阵K均值中心化处理 klX=kX-lmX*kX-kX*lmX+lmX*kX*lmX; %此函数是利用协方差矩阵对其进行主元分析。pk是主成分（特征向量构成的矩阵），lamd是核矩阵的协方差矩阵的特征值（特征值构成的向量。列向量），baifen是每个特征向量表征在观测量总方差中所占的百分数（向量。列向量）。 [pkX,lamdX,baifenX]=pcacov(cov(klX)); %对训练数据的特征向量进行归一化 for i=1:XX_row psX(:,i)=pkX(:,i)/sqrt(lamdX(i)); end %选择主元个数 num_pc=1;%主元个数的初值为1 while sum(lamdX(1:num_pc))/sum(lamdX)<0.85 num_pc=num_pc+1; end %求置信度分别在99%的T2统计量的T2控制限 t2a99=num_pc*(XX_row^2-1)*finv(0.99,num_pc,XX_row-num_pc)/(XX_row*(XX_row-num_pc)); % t2a95=num_pc*(XX_row-1)*icdf('f',0.95,num_pc,XX_row-num_pc)/(XX_row*(XX_row-num_pc)); %求置信度分别在99%的T2统计量的SPE控制限 tX=klX*psX;%计算SPE时要用到 for i=1:XX_row spe(i)=sum(tX(i,num_pc+1:XX_col).^2); end a=mean(spe); b=var(spe); Qa99=b/(2*a)*chi2inv(0.999,2*a^2/b); %% 在线监测 Xtest =dlmread('d01_te.dat');%载入检测数据 [XC_row,XC_col]=size(Xtest);%取测试数据阵的行和列数 % 用原始数据的均值和标准差对检测数据进行标准化 Xte_mean=repmat(Xt_mean,XC_row,1); Xte_std=repmat(Xt_std,XC_row,1); for i=1:XC_row for j=1:XX_col Xtest(i,j)=(Xtest(i,j)-Xte_mean(i,j))/Xte_std(i,j); end end %计算测试数据的核矩阵K for i=1:XC_row for j=1:XC_row kC(i,j)=exp(-(norm(Xtest(i,:)-Xtrain(j,:))).^2/400); end end lmC=1/XC_row*ones(XC_row);%定义测试数据的IN klC=kC-lmC*kX-kC*lmC+lmC*kX*lmC;%对测试数据的核矩阵K均值中心化处理 %T2统计量 tC=klC*psX; D=diag(lamdX); P=psX(:,1:num_pc); for i=1:XC_row T2(i)=klC(i,:)*P*inv(D(1:num_pc,1:num_pc))*P'*klC(i,:)'; end %SPE统计量 for i=1:XC_row Q(i)=sum(tC(i,num_pc+1:XC_col).^2); end %% 出图 %绘图 figure subplot(2,1,1); plot(1:XC_row,T2,'k'); title('对于故障1，用于检测的KPCA多变量统计量'); xlabel('采样数'); ylabel('T^2'); hold on; line([0,XC_row],[t2a99,t2a99],'LineStyle','--','Color','r'); subplot(2,1,2); plot(1:XC_row,Q,'k'); xlabel('采样数'); ylabel('SPE'); hold on; line([0,XC_row],[Qa99,Qa99],'LineStyle','--','Color','r');