failureprediction.m

这是一个TEP工业系统的仿真程序

%%------------------------PCA+CUSUM for failure prediction---------------------------------------------%%

close;
clc;
clear all;

%%----------------------get data for TE process------------------------------------------------------------------------------------------------------------------------%%
d00=load('d00.dat')';                          %% training file for the normal operating conditions, x=[XMEAS(1), XMEAS(2), ..., XMEAS(41), XMV(1), ..., XMV(11)]^T
% d00_te=load('d00_te.dat');                     %% testing file for the normal operating conditions
d03=load('d03.dat');                           %% training file for fault 3
% d03_te=load('d03_te.dat');                     %% testing file for fault 3
% d09=load('d09.dat');                           %% training file for fault 9
% d09_te=load('d09_te.dat');                     %% testing file for fault 9
%%----------------------------------------------------------------------------------------------------------------------------------------------------------------------%%

%%---------------------------PCA------22 XMEAS-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------%%
X=d00(:,1:22);
X=zscore(X);
[COEFF, SCORE, LATENT, TSQUARED] = princomp(X);   %%[1,2,3,4,5,6,9,10,11,13,14,16,17,18,21,22] LATENT is the eigenvalues of the covariance matrix of X, COEFF is the loadings, SCORE is the principal component scores, TSQUARED is the Hotelling's T-squared statistic for each observation in X.

%%-------------------choose k as the number of principal components----------------------------------------------------------------------------------------------------------------------------------%%
percent = input('确定方差贡献率限(0～1）：');                                                                                     %% alpha0 is the predetermined contribution rate, usually 85%
k=0;
for i=1:size(LATENT,1)                                                                            %% choose first k principal components
    
    alpha(i)=sum(LATENT(1:i))/sum(LATENT);
    
    if alpha(i)>=percent
        
        k=i;
        break;
        
    end
    
end
%%---------------------------------------------------------------------------------------------------------------------------------------------------------------%%

%%------------------------利用正态分布求SPE statistics---------------------------------------------------%%
Xp=zeros(size(X));             %% Xp is the estimation of original data using the PCA model, X=Xp+E

for j=1:k

    Xp=Xp+SCORE(:,j)*COEFF(:,j)';

end

for i=1:size(X,1)
    
    SPE(i)=sum((X(i,:)-Xp(i,:)).^2);
    
end

%% computing SPE limits
beta = input('确定统计阈值置信度(0～1）：'); %% beta is the inspection level
theta=zeros(3,1);
for i=1:3
    
    for j=k+1:size(X,2)
        theta(i)=theta(i)+LATENT(j)^(i);
    end

end

h0=1-2*theta(1)*theta(3)/3/(theta(2)^2);
SPEbeta=theta(1)*(norminv(beta)*(2*theta(2)*h0^2)^0.5/theta(1)+1+theta(2)*h0*(h0-1)/theta(1)^2)^(1/h0);

%%----------------------------------------------------------------------------------------------------%%

%%-----------------------------------利用F分布求T2 statistic--------------%%
T2knbeta=k*(size(X,1)-1)/(size(X,1)-k)*finv(beta,k,size(X,1));
%%-----------------------------------------------------------------------%%
%%--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------%%

%%---------------------T方统计量图----------------------------------------------%%
figure(1);
plot(TSQUARED);
hold on;
TS=T2knbeta*ones(size(X,1),1);
plot(TS,'r--');
title('T2 statistic of train data');
hold off;
%%-------------------------------------------------------------------%%

%%----------------------SPE统计量绘图---------------------------------------------%%
figure(2);
plot(SPE);
hold on;
S=SPEbeta*ones(size(X,1),1);
plot(S,'r--');
title('SPE statistic of train data');
hold off;
%%-------------------------------------------------------------------%%

%% %待测数据的预处理（步骤同样本数据）-------------------------------------------------------------------%%
XT=d03(:,1:22);
XT=zscore(XT);
XTb= XT*COEFF;

%% TS statistic figure
figure(3);
for i=1:size(XT,1);
    
    ts(i) = XTb(i,1:k) * diag( 1 ./ LATENT(1:k) ) * XTb(i,1:k)';

end
plot(ts);
title('Tsquare Control limits of test data');
hold on;
plot(T2knbeta(ones(1,size(X,1))),'r--');
hold off;

%% SPE statistic figure
figure(4);
for i=1:size(XT,1)
    
    spe(i) = XT(i,:)*(eye(size(XT,2))-COEFF(:,1:k)*COEFF(:,1:k)')*XT(i,:)';
        
end
plot(spe);
title('SPE Control limits of test data');
hold on;
plot(SPEbeta(ones(1,size(XT,1))),'r--');
hold off;
%%-------------------------------------------------------------------%%

%%-------------------------------------------------------------------%%
%%-------------------------------------------------------------------%%





























%%--------------------------------end-------------------------------------%%