基于kalman滤波的磨损预测算法matlab仿真

Posted 2022-11-18 fpga和matlab

目录

一、理论基础

二、核心程序

三、仿真测试结果

---

作者ID  ：fpga和matlab
CSDN主页：https://blog.csdn.net/ccsss22?type=blog
擅长技术：
1.无线基带,无线图传,编解码 
2.机器视觉,图像处理,三维重建 
3.人工智能,深度学习 
4.智能控制,智能优化
5.其他

一、理论基础

       磨损工程模型简介摩擦学的主要课题之一是对机器运行的可靠性和寿命作出评价,也就是对磨损作出预测,而磨损是一个复杂过程,它不仅与摩擦付的材料物理机械性能和化学性能有关,还与摩擦付的工况(压力,速度,温度),外部环境条件(周围介质,润滑剂,辐射等)和摩擦付结构参数有关,它在许多情况下包括化学浸蚀和物理损伤,工作条件的任何微小的改变都可能改变磨损过程的整个性质,由于这个理由,从实际"测试"中很难得到重复的或确实的结果,而甚至于对最周密地控制的实验室实验进行解释时也应特别。

      剩余寿命估计是工程系统预测与健康管理的关键． 目前，基于观测的系统退化数据进行剩余寿命估计得到了很大的关注． 由于系统随机退化过程和测量误差的影响，测量数据中不可避免包含退化随机性和测量不确定性． 然而，现有基于观测数据的剩余寿命估计研究中，没有将退化随机性和测量不确定性对估计的剩余寿命分布的影响同时考虑． 鉴于此，提出了一种基于 Wiener 过程且同时考虑随机退化和不确定测量的退化建模方法，利用 Kalman滤波技术，实现了潜在退化状态的实时估计． 在退化状态估计的基础上，得到了同时考虑退化状态不确定性和测量不确定性的解析剩余寿命分布． 此外，提出了一种基于极大似然方法的退化模型参数估计方法． 最后，通过陀螺仪的退化测量数据验证了本文提出的方法优于不考虑测量不确定性的方法，可以提高剩余寿估计的准确性．

        卡尔曼滤波是一种利用线性系统状态方程，通过系统输入输出观测数据，对系统状态进行最优估计的算法。由于观测数据中包括系统中的噪声和干扰的影响，所以最优估计也可看作是滤波过程。卡尔曼滤波（Kalman filtering）是一种利用线性系统状态方程，通过系统输入输出观测数据，对系统状态进行最优估计的算法。由于观测数据中包括系统中的噪声和干扰的影响，所以最优估计也可看作是滤波过程。
       数据滤波是去除噪声还原真实数据的一种数据处理技术，Kalman滤波在测量方差已知的情况下能够从一系列存在测量噪声的数据中，估计动态系统的状态。由于它便于计算机编程实现，并能够对现场采集的数据进行实时的更新和处理，Kalman滤波是目前应用最为广泛的滤波方法，在通信，导航，制导与控制等多领域得到了较好的应用。

 

 

二、核心程序

clc;
clear;
close all;
warning off;

rng(1);

................................................................
[R,C] = size(data);
Len   = 6;
%样本个数
Num   = R/Len;
Cnums = Len;
for i = 1:R/Len
    L1(:,i) = data(Len*(i-1)+1:Len*i,1); 
    L2(:,i) = data(Len*(i-1)+1:Len*i,3); 
    L3(:,i) = data(Len*(i-1)+1:Len*i,5); 
    L4(:,i) = data(Len*(i-1)+1:Len*i,7); 
    
    R1(:,i) = data(Len*(i-1)+1:Len*i,2); 
    R2(:,i) = data(Len*(i-1)+1:Len*i,4); 
    R3(:,i) = data(Len*(i-1)+1:Len*i,6); 
    R4(:,i) = data(Len*(i-1)+1:Len*i,8);     
end


figure;
subplot(221);
Y = L1;
plot(Y(1,:),'bo');
hold on
plot(Y(2,:),'bo');
hold on
plot(Y(3,:),'bo');
hold on
plot(Y(4,:),'bo');
hold on
plot(Y(5,:),'bo');
hold on
plot(Y(6,:),'bo');
hold on
grid on
axis([0,14,1,30]);
title('L1');

subplot(222);
Y = L2;
plot(Y(1,:),'bo');
hold on
plot(Y(2,:),'bo');
hold on
plot(Y(3,:),'bo');
hold on
plot(Y(4,:),'bo');
hold on
plot(Y(5,:),'bo');
hold on
plot(Y(6,:),'bo');
hold on
grid on
axis([0,14,1,30]);
title('L2');

subplot(223);
Y = L3;
plot(Y(1,:),'bo');
hold on
plot(Y(2,:),'bo');
hold on
plot(Y(3,:),'bo');
hold on
plot(Y(4,:),'bo');
hold on
plot(Y(5,:),'bo');
hold on
plot(Y(6,:),'bo');
hold on
grid on
axis([0,14,1,30]);
title('L3');

subplot(224);
Y = L4;
plot(Y(1,:),'bo');
hold on
plot(Y(2,:),'bo');
hold on
plot(Y(3,:),'bo');
hold on
plot(Y(4,:),'bo');
hold on
plot(Y(5,:),'bo');
hold on
plot(Y(6,:),'bo');
hold on
grid on
axis([0,14,1,30]);
title('L4');


figure;
subplot(221);
Y = R1;
plot(Y(1,:),'bo');
hold on
plot(Y(2,:),'bo');
hold on
plot(Y(3,:),'bo');
hold on
plot(Y(4,:),'bo');
hold on
plot(Y(5,:),'bo');
hold on
plot(Y(6,:),'bo');
hold on
grid on
axis([0,14,1,30]);
title('R1');

subplot(222);
Y = R2;
plot(Y(1,:),'bo');
hold on
plot(Y(2,:),'bo');
hold on
plot(Y(3,:),'bo');
hold on
plot(Y(4,:),'bo');
hold on
plot(Y(5,:),'bo');
hold on
plot(Y(6,:),'bo');
hold on
grid on
axis([0,14,1,30]);
title('R2');

subplot(223);
Y = R3;
plot(Y(1,:),'bo');
hold on
plot(Y(2,:),'bo');
hold on
plot(Y(3,:),'bo');
hold on
plot(Y(4,:),'bo');
hold on
plot(Y(5,:),'bo');
hold on
plot(Y(6,:),'bo');
hold on
grid on
axis([0,14,1,30]);
title('R3');

subplot(224);
Y = R4;
plot(Y(1,:),'bo');
hold on
plot(Y(2,:),'bo');
hold on
plot(Y(3,:),'bo');
hold on
plot(Y(4,:),'bo');
hold on
plot(Y(5,:),'bo');
hold on
plot(Y(6,:),'bo');
hold on
grid on
axis([0,14,1,30]);
title('R4');


%%
%去除异常点后的L1均值数据,注意，你里面有一堆错误采集数据，需要根据实际情况来剔除，然后再进行后续的计算，所以下面这个人工剔除的过程，需要自己修改
Y = L1;
for i = 1:13
    datas(i) = (Y(1,i) + Y(2,i) + Y(3,i) + Y(4,i) + Y(5,i) + Y(6,i))/6;
    if i == 4
       datas(i) = (Y(1,i) + Y(2,i) + Y(3,i) + Y(4,i) + Y(5,i))/5; 
    end
    if i == 7
       datas(i) = (Y(1,i) + Y(2,i) + Y(3,i) + Y(4,i) + Y(6,i))/5;  
    end    
    if i == 12
       datas(i) = (Y(1,i) + Y(2,i) + Y(3,i) + Y(6,i))/4; 
    end
    if i == 13
       datas(i) = (Y(2,i) + Y(3,i) + Y(5,i) + Y(6,i))/4;  
    end
end




figure;
plot(datas,'r-o');
hold on
Y = L1;
plot(Y(1,:),'bo');
hold on
plot(Y(2,:),'bo');
hold on
plot(Y(3,:),'bo');
hold on
plot(Y(4,:),'bo');
hold on
plot(Y(5,:),'bo');
hold on
plot(Y(6,:),'bo');
hold on
grid on
axis([0,14,1,30]);
title('L1');
xlabel('Times');
grid on

%%
data_out2 = func_kalman_predict(datas,25);
figure;
subplot(121);
plot(data_out2,'r-s');
hold on
plot(datas,'b-*');
legend('预测值','原始测量值均值（剔除异常点）');
xlabel('月份');
ylabel('thickness/mm');
grid on
axis([0,30,0,20]);

subplot(122);
plot(data_out2,'r-s');
hold on
Y = L1;
plot(Y(1,:),'bo');
hold on
plot(Y(2,:),'bo');
hold on
plot(Y(3,:),'bo');
hold on
plot(Y(4,:),'bo');
hold on
plot(Y(5,:),'bo');
hold on
plot(Y(6,:),'bo');
hold on
legend('预测值','原始测量值');
xlabel('月份');
ylabel('thickness/mm');
grid on
axis([0,30,0,20]);

figure;
plot((data_out2(1:length(datas))-datas).^2,'b-*');
grid on
xlabel('月份');
ylabel('均方误差/mm');



%计算剩余寿命
for i = 10:100
    data_out2 = func_kalman_predict(datas,i);
    if data_out2(end) < 1
       LIFE = i;
       break
    end
end
disp('总寿命月份');
LIFE
figure;
plot(data_out2,'r-s');
hold on
plot(datas,'b-*');
legend('预测值','原始测量值均值（剔除异常点）');
xlabel('月份');
ylabel('thickness/mm');
grid on
axis([0,LIFE,0,20]);
title(['总寿命：',num2str(LIFE)]);


for i = 1:length(data_out2)-1
    ms(i) = abs(data_out2(i+1)-data_out2(i));
end
%损耗期望
E=mean(ms)
%损耗均方误差
D=std(ms)

 

三、仿真测试结果

 

 

 A28-64