close all;
clear all;
clc;
%% 正常数据
Xtrain =dlmread('d00_te.dat');
X_mean = mean(Xtrain); %平均值
X_std = std(Xtrain); %方差
[X_row ,X_col]= size(Xtrain); %训练集矩阵的行数和列数
for i = 1:X_col
Xtrain(:,i) = (Xtrain(:,i)-X_mean(i))./X_std(i); %训练集矩阵每一列减去平均值除以方差
end
[U,S,V]=svd(Xtrain./sqrt(size(Xtrain,1)-1));
D= S.^2;
lamda=diag(D);
num_pc=1;
while sum(lamda(1:num_pc))/sum(lamda)<0.85
num_pc=num_pc+1;
end
P=V(:,1:num_pc);%取与lamda相对应的特征向量
%求置信度为`99%时的T2统计控制限
T2UCL1=num_pc*(X_row-1)*(X_row+1)*finv(0.99,num_pc,X_row - num_pc)/(X_row*(X_row - num_pc));
%置信度为99%的Q统计控制限
for i = 1:3
theta(i) = sum((lamda(1:num_pc)).^2i);
end
h0 = 1 - 2*theta(1)*theta(3)/(3*theta(2)^2);
ca = norminv(0.99,0,1);
QUCL = theta(1)*(h0*ca*sqrt(2*theta(2))/theta(1) + 1 + theta(2)*h0*(h0 - 1)/theta(1)^2)^(1/h0);
%% 在线监测
Xtest =dlmread('d01_te.dat');
Xte_mean=repmat(X_mean,size(Xtest,1),1);
Xte_std=repmat(X_std,size(Xtest,1),1);
for i=1:size(Xtest,1)
for j=1:size(Xtest,2)
Xtest(i,j)=(Xtest(i,j)-Xte_mean(i,j))/Xte_std(i,j);
end
end
[X_row,X_col]= size(Xtest);
lamda=diag(lamda);
D=lamda(1:num_pc,1:num_pc);
%计算T2统计量
T2 = zeros(X_row,1);
for i=1:X_row
T2(i)=Xtest(i,:)*P*inv(D)*P'*Xtest(i,:)';
end
%计算SPE统计量
SPE= zeros(X_row,1);
P=V(:,num_pc+1:X_col);
[r,y]=size(P*P');
I=eye(r,y);
for l=1:X_row
SPE(l)= Xtest(l,:)*(I - P*P')*(I - P*P')'*Xtest(l,:)';
end
%%
%绘图
figure
subplot(2,1,1);
plot(1:X_row,T2,'k');
title('对于故障1,用于检测的PCA多变量统计量');
xlabel('采样数');
ylabel('T^2');
hold on;
line([0,X_row],[T2UCL1,T2UCL1],'LineStyle','--','Color','r');
subplot(2,1,2);
plot(1:X_row,SPE,'k');
xlabel('采样数');
ylabel('SPE');
hold on;
line([0,X_row],[QUCL,QUCL],'LineStyle','--','Color','r');