test_TE.asv
上传用户:nxf103
上传日期:2017-12-04
资源大小:2690k
文件大小:3k
开发平台:

Matlab

  1. global ks;%核函数半径
  2. load x00;
  3. load t02;
  4. ks=50;
  5. train_num = 960;% 训练样本数
  6. fault_num = 800;
  7. variable_num=16;
  8. PCA_num=5;%选取最大主元数
  9. K=2;%KPCA选取最大主元数
  10. %%%%%对建模样本作归一化处理%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  11. [n,m]=size(x00);
  12. E=mean(x00);
  13. Train=zeros(n,m);
  14. for j=1:m,
  15.    V(j)=0; 
  16.    for i=1:n,
  17.      V(j)=V(j)+(x00(i,j)-E(j))^2;
  18.    end
  19.    V(j)=sqrt(V(j)/(n-1));
  20.    for i=1:n,
  21.        Train(i,j)=(x00(i,j)-E(j))/V(j);
  22.    end 
  23. end
  24. %%%%%对故障测试样本作归一化处理%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  25.  [n1,m1]=size(t02);
  26.  fault=zeros(n1,m1); 
  27.  for j=1:m1,
  28.     for i=1:n1,
  29.         fault(i,j)=(t02(i,j)-E(j))/V(j);
  30.     end 
  31.  end
  32.  %%%%% %主元分析%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  33.   [pc,score,latent,tsquare]=princomp(Train);
  34.   PCA=pc(:,1:PCA_num);
  35.   score_PCA=score(:,1:PCA_num);
  36.   latent_PCA=diag(latent(1:PCA_num));
  37. %%%%% %新测试数据(故障数据)作主元分析%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  38.   PCA_Test=zeros(n1,PCA_num);
  39.   Fault_guji=zeros(n1,m1);
  40.   TT=zeros(1,n1);
  41.   
  42. for i=1:n1,
  43.      for j=1:PCA_num,
  44.        PCA_Test(i,j)=fault(i,:)*PCA(:,j); 
  45.        Fault_guji(i,:)=Fault_guji(i,:)+fault(i,:)*PCA(:,j)*PCA(:,j)';
  46.      end
  47.   %求T2统计
  48.   TT(i)=PCA_Test(i,:)*inv(latent_PCA)*PCA_Test(i,:)';
  49.   %求SPE统计
  50. end 
  51.  Fault_ee=fault-Fault_guji;
  52.   for i=1:n1,
  53.   Fault_SPE(i)=Fault_ee(i,:)*Fault_ee(i,:)';
  54.   end
  55.    figure(1);  
  56.    subplot(2,1,1);
  57.    plot(1:960,TT,1:8:960,20,'k--');
  58.    xlabel('Sample Number');
  59.    ylabel('PCA-T2');
  60.    subplot(2,1,2);
  61.    plot(1:960,Fault_SPE,1:8:960,30,'k--');
  62.    xlabel('Sample Number');
  63.    ylabel('PCA-SPE'); 
  64. %%%%% %Kernel主元分析%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
  65. x00=[x00(1,:); x00(80,:);x00(120,:)];
  66. [tr,meanp,stdp] = prestd(x00');
  67. te = trastd(t02',meanp,stdp);
  68. tr=tr';
  69. te=te';
  70. trd=data(tr);
  71. ted=data(te);
  72. %[r,a]=train(kpca,tr);
  73. [r,a]=train(kpca({kernel('rbf',ks)}),trd);
  74. d=test(a,ted);
  75. val=a.e_val;
  76. vec=a.e_vec;
  77. dvec=d.X;
  78. s=inv(diag(val(1:K)));
  79. for i=1:960
  80.   tsquare(i)=dvec(i,1:K)*s*dvec(i,1:K)';
  81.   SPE(i)=dvec(i,:)*dvec(i,:)'-dvec(i,1:K)*dvec(i,1:K)';
  82. end
  83. figure(2);
  84. subplot(2,1,1);
  85. plot(1:960,tsquare,1:8:960,0.025,'k--');
  86. xlabel('Sample Number(N)');
  87. ylabel('KPCA-T2');
  88. subplot(2,1,2);
  89. plot(1:960,SPE,1:8:960,0.025,'k--');
  90. xlabel('Sample Number(N)');
  91. ylabel('KPCA-SPE');
  92. %[cn,CSPE]=kpcak(te,tr,vec(:,1:K),s,ks);
  93. %figure(3)
  94. subplot(2,1,1);
  95. plot(cn);
  96. xlabel('Sample Number(N)');
  97. ylabel('CT2');
  98. subplot(2,1,2);
  99. plot(CSPE);
  100. xlabel('Sample Number(N)');
  101. ylabel('CSPE');