文章标题

MSET的MATLAB实现

具体流程如下所示：分别导入历史状态集合和观测状态集合——数据归一化——求出预测值——作图，画出各个变量的观测值和预测值——用欧氏距离表达他们的距离——滑动窗口统计法求相似度曲线和监测线。

clear
clc
[a1]=xlsread('shuju','F3:F1322');%发电功率
[b1]=xlsread('shuju','I3:I1322');%2号引风机入口烟气压力
[c1]=xlsread('shuju','O3:O1322');%2号引风机轴承1振动
[d1]=xlsread('shuju','L3:L1322');%2号引风机轴承1温度
NUM=length(a1);
%归一化处理
a=zeros(1,NUM);
b=zeros(1,NUM);
c=zeros(1,NUM);
d=zeros(1,NUM);
a1=a1';
b1=b1';
c1=c1';
d1=d1';
for m=1:NUM
    a(1,m)=(a1(1,m)-min(a1))/(max(a1)-min(a1));
    b(1,m)=(b1(1,m)-min(b1))/(max(b1)-min(b1));
    c(1,m)=(c1(1,m)-min(c1))/(max(c1)-min(c1));
    d(1,m)=(d1(1,m)-min(d1))/(max(d1)-min(d1));
end
D=[a;b;c;d]; %归一化后历史向量集合
D_=[a1;b1;c1;d1];%真实的历史向量集合
[x]=xlsread('shuju','F1331:F1470');
[y]=xlsread('shuju','I1331:I1470');
[z]=xlsread('shuju','O1331:O1470');
[p]=xlsread('shuju','L1331:L1470');
num=length(x);
%将读取的数据倒置
x=x';
y=y';
z=z';
p=p';
%归一化处理
x1=zeros(1,num);
y1=zeros(1,num);
z1=zeros(1,num);
p1=zeros(1,num);
for m=1:num
    x1(1,m)=(x(1,m)-min(x))/(max(x)-min(x));
    y1(1,m)=(y(1,m)-min(y))/(max(y)-min(y));
    z1(1,m)=(z(1,m)-min(z))/(max(z)-min(z));
    p1(1,m)=(p(1,m)-min(p))/(max(p)-min(p));
end
Kobs=[x1;y1;z1;p1];%归一化后的观测值
Kobs_=[x;y;z;p];%真实的观测值
%MEST模型
Temp=zeros(NUM);
for i=1:NUM
    for j=1:NUM
         Temp(i,j)=norm(D_(:,i)-D_(:,j),2); 
    end
end
Temp1=ones(NUM,num);
for ii=1:NUM
    for jj=1:num
         Temp1(ii,jj)=norm(D_(:,ii)-Kobs_(:,jj),2); 
    end
end
Kest_=D_*inv(Temp)*Temp1;%真实的预测值
%预测值归一化
Kest=zeros(4,num);
 for k=1:num
    Kest(1,k)=(Kest_(1,k)-min(Kest_(1,:)))/(max(Kest_(1,:))-min(Kest_(1,:)));
    Kest(2,k)=(Kest_(2,k)-min(Kest_(2,:)))/(max(Kest_(2,:))-min(Kest_(2,:)));
    Kest(3,k)=(Kest_(3,k)-min(Kest_(3,:)))/(max(Kest_(3,:))-min(Kest_(3,:)));
    Kest(4,k)=(Kest_(4,k)-min(Kest_(4,:)))/(max(Kest_(4,:))-min(Kest_(4,:)));
 end
err=Kobs_-Kest_;
xerror=err(1,:);
yerror=err(2,:);
zerror=err(3,:);
perror=err(4,:);
%画图：
outputpict(Kest_,Kobs_,xerror,yerror,zerror,perror);

%求相似度曲线--欧氏距离
sim=zeros(1,num);
for i=1:num
    sim(1,i)=1-norm(Kobs(:,i)-Kest(:,i),2)/0.75;
end
% figure(5)
% plot(sim,'LineWidth',2);
%求相似度检测线--滑动窗口统计法
N=30;%N为窗口宽度
sim_js=zeros(1,num+1-N);
sim_jssum=zeros(1,num+1-N);
for i=1:num+1-N
    for j=0:N-1
        sim_jssum(1,i)=sim_jssum(1,i)+sim(i+j);
    end
    sim_js(1,i)=sim_jssum(1,i)/N;
end
sim_jsmin=min(sim_js);
sen=0.98;%sen为预警的灵敏性
figure(6)
plot(sim_js,'LineWidth',2);
hold on
plot([0,num],[0.66869,0.66869],'r','LineWidth',2);%从健康数据验证得到检测线为0.66869

数据用的是引风机的数据，监测线由健康数据求出，通过模拟温度的故障可以检验模型是否能发出预警。其中作图的函数如下：

function outputpict(Kest_,Kobs_,xerror,yerror,zerror,perror)
figure(1)
subplot(2,1,1)
plot(Kest_(1,:),'k','LineWidth',1.5);
hold on
plot(Kobs_(1,:),'--*r','LineWidth',1.2);
title('发电功率的预测值与观测值的比较');
grid on;
subplot(2,1,2)
plot(xerror,'LineWidth',1.5);
title('发电功率的预测值与观测值之间的误差');
grid on;
%输出y
figure(2)
subplot(2,1,1)
plot(Kest_(2,:),'k','LineWidth',1.5);
hold on
plot(Kobs_(2,:),'--*b','LineWidth',1.2);
title('入口烟气压力的预测值与观测值的比较');
grid on;
subplot(2,1,2)
plot(yerror,'LineWidth',1.5);
grid on;
title('入口烟气压力的预测值与观测值之间的误差');
%输出z
figure(3)
subplot(2,1,1)
plot(Kest_(3,:),'k','LineWidth',1.5);
hold on
plot(Kobs_(3,:),'--*g','LineWidth',1.2);
title('引风机轴承1水平振动的预测值与观测值的比较');
grid on;
subplot(2,1,2)
plot(zerror,'LineWidth',1.5);
title('引风机轴承1水平振动的预测值与观测值之间的误差');
grid on;
%输出p
figure(4)
subplot(2,1,1)
plot(Kest_(4,:),'k','LineWidth',1.5);
hold on
plot(Kobs_(4,:),'--*m','LineWidth',1.2);
title('引风机轴承1温度的预测值与观测值的比较');
grid on;
subplot(2,1,2)
plot(perror,'LineWidth',1.5);
title('引风机轴承1温度的预测值与观测值之间的误差');
grid on;
end