PCA故障诊断matlab实现200654

matlab 实现PCA故障诊断功能，带测试测试数据，可直接运行

%建立模型：
Xtrain = CC（1:1000:）;
Xtrain_Store = Xtrain;
%测试数据
Xtest = CC（1001:1500:）;
Xtest_Store = Xtest;
size_= size（Xtest）;
test_Count = size_（1）;
%总输入变量个数
in_count = 10;
%标准化处理：
X_mean = mean（Xtrain）;  %按列求Xtrain平均值                           
X_std = std（Xtrain）;    %求标准差                      
[X_rowX_col] = size（Xtrain）; %求Xtrain行、列数               
% for i = 1:X_col
%Xtrain（:i）=（Xtrain（:i） - X_mean（i）./X_std（i））;
%Xtest（:i） = （Xtest（:i） - X_mean（i）./X_std（i））;
% end                                                
Xtrain=（Xtrain-repmat（X_meanX_row1））./repmat（X_stdX_row1）;

%求协方差矩阵
sigmaXtrain = cov（Xtrain）;
%对协方差矩阵进行特征分解，lamda为特征值构成的对角阵，T的列为单位特征向量，且与lamda中的特征值一一对应：
[Tlamda] = eig（sigmaXtrain）;                            
% disp（‘特征根（由小到大）‘）;
% disp（lamda）;
% disp（‘特征向量：‘）;
% disp（T）;                                            

%取对角元素（结果为一列向量），即lamda值，并上下反转使其从大到小排列，主元个数初值为1，若累计贡献率小于90%则增加主元个数
D = flipud（diag（lamda））;                            
num_pc = 1;                                         
while sum（D（1:num_pc））/sum（D） < 0.9   
num_pc = num_pc +1;
end                                                 

%取与lamda相对应的特征向量
P = T（:X_col-num_pc+1:X_col）;                            

%求置信度为99%、95%时的T2统计控制限                       
T2UCL1=num_pc*（X_row-1）*（X_row+1）*finv（0.999num_pcX_row - num_pc）/（X_row*（X_row - num_pc））;
T2UCL2=num_pc*（X_row-1）*（X_row+1）*finv（0.95num_pcX_row - num_pc）/（X_row*（X_row - num_pc））;

%置信度为99%的Q统计控制限
for i = 1:3
    theta（i） = sum（（D（num_pc+1:X_col））.^i）;
end
h0 = 1 - 2*theta（1）*theta（3）/（3*theta（2）^2）;
ca = norminv（0.99901）; %正态累积分布的反函数
QUCL = theta（1）*（h0*ca*sqrt（2*theta（2））/theta（1） + 1 + theta（2）*h0*（h0 - 1）/theta（1）^2）^（1/h0）;                           

%测试：
%标准化处理
n = size（Xtest1）;
Xtest=（Xtest-repmat（X_meann1））./repmat（X_stdn1）;

%求T2统计量，Q统计量
[ry] = size（P*P‘）;
I = eye（ry）;

T2 = zeros（n1）;
Q = zeros（n1）;
for i = 1:n
    T2（i）=Xtest（i:）*P*inv（lamda（in_count-num_pc+1:in_countin_count-num_pc+1:in_count））*P‘*Xtest（i:）‘;                                           
    Q（i） = Xtest（i:）*（I - P*P‘）*Xtest（i:）‘;                                                                                    
end

%绘图
    %figure
    subplot（211）;
    plot（1:nT2‘k‘）;                                    
    title（‘主元分析统计量变化图‘）;
    xlabel（‘采样数‘）;
    ylabel（‘T^2‘）;
    hold on;
    line（[0n][T2UCL1T2UCL1]‘Linestyle‘‘--‘‘Color‘‘r‘）;
    line（[0n][T2UCL2T2UCL2]‘Linestyle‘‘--‘‘Color‘‘g‘）;
    
    subplot（212）;
    plot（1:nQ‘k‘）;
    xlabel（‘采样数‘）;
    ylabel（‘Q‘）;
    hold on;
    line（[0n][QUCLQUCL]‘Linestyle‘‘--‘‘Color‘‘r‘）;
    
%贡献图
%1.确定造成失控状态的得分
S = Xtest（test_Count:）*P（:1:num_pc）;
r = [ ];
for i = 1:num_pc
    if S（i）^2/lamda（i） > T2UCL1/num_pc
        r = cat（2ri）;
    end
end

%2.计算每个变量相对于上述失控得分的贡献
cont = zeros（length（r）in_count）;
for i = length（r）
    for j = 1:in_count

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3756  2014-04-25 16:32  Test2.m

     文件     355986  2014-04-25 16:34  PCAdata.mat

----------- ---------  ---------- -----  ----

               359742                    2