故障分析基于matlab ICA故障监测含Matlab源码 1590期

Posted 2022-01-03 紫极神光

一、获取代码方式

获取代码方式1：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

获取代码方式2：
完整代码已上传我的资源：【故障分析】基于matlab ICA故障监测【含Matlab源码 1590期】

备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、部分源代码

clear all
clc
cont=0.9;
alpha=0.99;
c_aphla=2.32;    %对应于0.99，置信水平为95%对应分位点为1.645
%建模
xn=textread('d10.dat');
x=zscore(xn);
x=x';
[m,n]=size(x);
[S,Q,B,evals,evecs]=ICA_normal(x)
[ICn,Bn,d,ICa]=sort_IC(S,Q,B,cont);
[I2,SPE]=variable_c(x,ICn,Q,Bn);
[f1,x1,u1]=ksdensity(I2);%I2单变量核密度估计
ConInt1=ComCon(f1,x1,alpha);
I2_limit=ConInt1(2);
SPE_limit=ksdensity(SPE,alpha,'function','icdf');

%在线监控
Xn=textread('d10_te.dat');
X=zscore(Xn);
X=X';
fault_I2_num=[];
fault_SPE_num=[];
[S_new,Q_new,B_new]=ICA_monitor(X,evals,evecs);
[ICn_new,Bn_new,d_new,ICa_new]=sort_IC(S_new,Q_new,B_new,cont);
[I2_new,SPE_new]=variable_c(X,ICn_new,Q_new,Bn_new);
for i=1:n
    if I2_new(i)>I2_limit
        fault_I2_num=[fault_I2_num,i];
    end;
    if SPE_new(i)>SPE_limit
        fault_SPE_num=[fault_SPE_num,i];
    end;
end;
fault_I2_num
fault_SPE_num

figure(1)
plot(1:length(I2_new),I2_new,'b');
hold on
plot(1:length(I2_new),ones(length(I2_new),1)*I2_limit,'r-');
xlabel('样本');ylabel('样本点的I2贡献值');
legend('I2统计量贡献','I2统计量控制限');
figure(2)
plot(1:length(SPE_new),SPE_new,'b');
hold on
plot(1:length(SPE_new),ones(length(SPE_new),1)*SPE_limit,'r-');
legend('SPE统计量贡献','SPE统计量控制限');
xlabel('样本');ylabel('样本点的SPE贡献值');

%故障诊断
new1=fault_I2_num(1);
new2=fault_SPE_num(1);
Xc_new=inv(Q_new)*Bn_new*ICn_new;
cont_I2=Xc_new(:,new1)*norm(S_new(:,new1))/norm(Xc_new(:,new1));
cont_SPE=(X(:,new2)-Xc_new(:,new2)).^2;
figure(3)
bar(cont_I2);
xlabel('变量');ylabel('变量对I2贡献');
figure(4)
bar(cont_SPE);
xlabel('变量');ylabel('变量对SPE贡献');
%对已知的密度函数求取控制限
function ConInt=ComCon(density,xmesh,alpha)	%density 密度矩阵 
							%xmesh为密度对应的x值的行向量（或矩阵）
							%alpha 置信度
ConInt=zeros(1,2);			%输出ConInt为控制限
n=length(density);			%计算每一列的控制限

x1=min(xmesh);
x2=max(xmesh);
dx=(x2-x1)/(n-1);

%黄金分割搜索从开始到某一点的面积为(1-alpha)/2,以求取下分位数;
while(1)
     x3=x1+(sqrt(5)-1)/2*(x2-x1);
     t=x1:dx:x3;
%    s=trapz(t,density(1:length(t)));
     s=trapz(density(1:length(t)))*dx;
     if abs(x2-x1)<1.0e-3 
%       if s-alpha<1.0e-4
        ConInt(1)=x3; 
     break;
%    elseif s<alpha
     elseif s<(1-alpha)/2
            x1=x3;
     else
          x2=x3;
     end
end
%     求上分为数，从开始到某一点的面积为1-（1-alpha)/2;
x1=min(xmesh);
x2=max(xmesh);
while(1)
      x3=x1+(sqrt(5)-1)/2*(x2-x1);
      t=x1:dx:x3;
%     s=trapz(t,density(1:length(t)));
      s=trapz(density(1:length(t)))*dx;
      if abs(x2-x1)<1.0e-3 
         ConInt(2)=x3; 
         break;
      elseif s<(1+alpha)/2
             x1=x3;
      else
           x2=x3;
      end
end

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 沈再阳.精通MATLAB信号处理[M].清华大学出版社，2015.
[2]高宝建,彭进业,王琳,潘建寿.信号与系统——使用MATLAB分析与实现[M].清华大学出版社，2020.
[3]王文光,魏少明,任欣.信号处理与系统分析的MATLAB实现[M].电子工业出版社，2018.
[4]严华,申雨.基于Matlab的轴承故障诊断分析[J].中国水运(下半月). 2021,21(02)