normalizedmissfault.m

基于粒子滤波的故障检测

close all;
clear all;
n = 1:600;%sample steps
stdw = sqrt(10);
ngrid = 50;
npar = 500;%particle number
N = length(n);

zideal=zeros(1,N);
error=zeros(1,N);
d=zeros(1,N);
Dfault=zeros(1,N);
M=20;

H=0.1:0.1:0.9;
faultalarm=zeros(1,9);

A=50;F=300;D=300;

for h=1:9  %threshold from 50,60,70,......200
  
   faultN=0;
   for s=1:50
%N=500
% Generate the state process and observations
x0 = 0.1;
c0 = 1;
b0=25;
xpath = zeros(1,N);
xmean=0;xvariance=0.1;
ymean=0;yvariance=1;
xnoise=gauss(xmean,xvariance,N); 
ynoise=gauss(ymean,yvariance,N); 

b=b0;
xpath(1) = x0/2 + b*x0/(1+x0^2) +8*cos(0) +xnoise(1) ;
for i=2:N
    if i<301
        b=b0;
    else 
        b=b0*5;
    end
xpath(i) = xpath(i-1)/2 + b*xpath(i-1)/(1+xpath(i-1)^2) +8*cos(1.2*(i-1)) + xnoise(i);
end
zpath = 1/20*(xpath.^2) + ynoise;
% Particle filter with resampling
w = ones(npar,1)/npar;
xprev = randn(npar, 1);
SParMat = zeros(npar, N);
SXParMat = zeros(npar, N);
sxparpath = zeros(1,N);


for i=1:N
    faultnumber=0;
xnext = drawpar(xprev, stdw, i);
xs = (xnext.^2)/20;
w = w.*(1/sqrt(2*pi)*exp(-((zpath(i)-xs).^2)/2));
l=1/sqrt(2*pi)*exp(-((zpath(i)-xs).^2)/2);
%Li=sum(l)/npar;%i时刻的所有粒子似然函数值均值
%likelihood(i)=Li;
%Di=0;
%if i<20
 %   for j=1:i
  %      Di=Di+(-(log(likelihood(j))));
   % end
%else 
 %   for j=i-20+1:i
  %      Di=Di+(-(log(likelihood(j))));
   % end
%end
%D(i)=Di;
w = w/sum(w);
SParMat(:,i) = w;
SXParMat(:,i) = xnext;
sxparpath(i) = w'*xnext;
[xprev, w] = impResample(xnext, w);
zideal(i)=1/20*sxparpath(i)^2;%the ideal observation at time i
error(i)=abs(zideal(i)-zpath(i));
Error=0;
if i<M
    for j=1:i
        Error=Error+error(j);
    end
    d(i)=1/i*Error;
else
    for j=i-M+1:i
        Error=Error+error(j);
    end
    d(i)=1/M*Error;
end
%if i>300
 %   if d(i)<h
  %      missnumber=missnumber+1;
   % else
    %    missnumber=missnumber;
    %end
%end
end
a=sum(d);
Dfault(i)=d(i)/a;
for i=1:300
    if Dfault(i)>h*0.1
        faultnumber=faultnumber+1;
    else
        faultnumber=faultnumber;
    end
end
faultN=faultN+faultnumber;
   end
faultalarm(h)=faultN/(A*F);
end




%normalized missalarm 
%to produce the missing alarm of SIR state estimation and smoothed residual

n = 1:600;%sample steps
stdw = sqrt(10);
ngrid = 50;
npar = 500;%particle number
N = length(n);

zideal=zeros(1,N);
error=zeros(1,N);
d=zeros(1,N);
Dmiss=zeros(1,N);
M=20;

H=0.1:0.1:0.9;
missalarm=zeros(1,9);
A=50;F=300;D=300;

for h=1:9  %threshold from 50,60,70,......200
  
   missN=0;
   for s=1:50
%N=500
% Generate the state process and observations
x0 = 0.1;
c0 = 1;
b0=25;
xpath = zeros(1,N);
xmean=0;xvariance=0.1;
ymean=0;yvariance=1;
xnoise=gauss(xmean,xvariance,N); 
ynoise=gauss(ymean,yvariance,N); 

b=b0;
xpath(1) = x0/2 + b*x0/(1+x0^2) +8*cos(0) +xnoise(1) ;
for i=2:N
    if i<301
        b=b0;
    else 
        b=b0*5;
    end
xpath(i) = xpath(i-1)/2 + b*xpath(i-1)/(1+xpath(i-1)^2) +8*cos(1.2*(i-1)) + xnoise(i);
end
zpath = 1/20*(xpath.^2) + ynoise;
% Particle filter with resampling
w = ones(npar,1)/npar;
xprev = randn(npar, 1);
SParMat = zeros(npar, N);
SXParMat = zeros(npar, N);
sxparpath = zeros(1,N);


for i=1:N
    missnumber=0;
xnext = drawpar(xprev, stdw, i);
xs = (xnext.^2)/20;
w = w.*(1/sqrt(2*pi)*exp(-((zpath(i)-xs).^2)/2));
l=1/sqrt(2*pi)*exp(-((zpath(i)-xs).^2)/2);
%Li=sum(l)/npar;%i时刻的所有粒子似然函数值均值
%likelihood(i)=Li;
%Di=0;
%if i<20
 %   for j=1:i
  %      Di=Di+(-(log(likelihood(j))));
   % end
%else 
 %   for j=i-20+1:i
  %      Di=Di+(-(log(likelihood(j))));
   % end
%end
%D(i)=Di;
w = w/sum(w);
SParMat(:,i) = w;
SXParMat(:,i) = xnext;
sxparpath(i) = w'*xnext;
[xprev, w] = impResample(xnext, w);
zideal(i)=1/20*sxparpath(i)^2;%the ideal observation at time i
error(i)=abs(zideal(i)-zpath(i));
Error=0;
if i<M
    for j=1:i
        Error=Error+error(j);
    end
    d(i)=1/i*Error;
else
    for j=i-M+1:i
        Error=Error+error(j);
    end
    d(i)=1/M*Error;
end
%if i>300
 %   if d(i)<h
  %      missnumber=missnumber+1;
   % else
    %    missnumber=missnumber;
    %end
%end
end
D=sum(d);
Dmiss(i)=d(i)/D;
for i=301:N
    if Dmiss(i)<h*0.1
        missnumber=missnumber+1;
    else
        missnumber=missnumber;
    end
end
missN=missN+missnumber;
   end
missalarm(h)=missN/(A*F);
end



figure;
plot(H,faultalarm,'g-*');hold on;
plot(H,missalarm,'r-o');
title('fault alarm and miss alarm based on state estimation and smoothed residual');
legend('faultalarm','missalarm');