链接:https://pan.baidu.com/s/1s6T7UpY1G2m-A1876etmBg
提取码:1234
一、参考资料
1、基于强跟踪滤波器的在线故障诊断方法
2、线性连续系统的离散化
二、理论推导
画重点:当考虑系统受到载荷激励时,微分方程右边为载荷激励;当考虑系统受到位移和速度的激励时,右边为位移速度激励以刚度k和阻尼c为系数的线性合作。这是两个不同的动力学系统。如果要卡尔曼滤波监测k和c,则需要扩展k和c两个状态变量,这时只能用第二种系统模型,否则卡尔曼滤波在状态预测时,不准确的k和c的估计值会影响到载荷激励,与真实的载荷差别很大,使估计误差的方差不仅不减小,反而会增大,致使算法无法收敛。模型错误是卡尔曼滤波不收敛的重要原因。
三、线性系统离散及验证
四、强跟踪卡尔曼滤波故障诊断
假设系统初始刚度为7.55/(pi/180),0.3秒后系统发生故障,刚度变成450/(pi/180),刚度初始估计值为0.5*7.55/(pi/180),看卡尔曼滤波可否快速跟踪刚度的变化从而进行故障判断。
隔振系统仿真代码
%kalman_sim.m
%Plant
function [u]=kalman_sim(u1,u2,u3,u4,u5)
%u1:输入飞轮位移
s_flywheel=u1;
%u2:输入飞轮转速
v_flywheel=u2;
%u3:系统噪声W
W=u3;
%u4:测量噪声V
V=u4;
%u5:时间t
t=u5;
%s一轴转角 v一轴转速 m一轴当量质量 k离合器刚度 c离合器阻尼
%Ts采样时间(由输入信号采样频率确定)
%F_s_flywheel飞轮转角 F_v_flywheel飞轮转速
persistent s v m k c Ts
if t==0
fs=1792.1;
Ts=1/fs;
s=0.0;
v=0.0;
m=0.3511;%kgm^2
k=7.55/(pi/180);%Nm/rad
wn=sqrt(k/m);
s_=0.25;
%wd=sqrt(1-s_^2)*wn;
c=2*s_*wn*m;
z=v+V;
end
if t>0.3
k=450/(pi/180);%假设5秒后刚度出现突变为450
end
if t>0
F_input=c*v_flywheel+k*s_flywheel;
s=s+Ts*v;
v=-k/m*Ts*s+(1-c/m*Ts)*v+Ts/m*F_input+Ts*W;
z=v+V;
end
u(1)=v;
u(2)=k;
u(3)=c;
u(4)=z;
强跟踪卡尔曼滤波代码
%kalman_k_c.m
%Plant
function [u]=kalman_sim(u1,u2,u3,u4)
%u1:输入飞轮位移
s_flywheel=u1;
%u2:输入飞轮转速
v_flywheel=u2;
%u3:系统噪声W
z=u3;
%u4:时间t
t=u4;
%persistent v m i
persistent x Ts m Q R P H Jf V0
if t==0
fs=1792.1;
Ts=1/fs;
m=0.3511;%kgm^2
k=7.55/(pi/180);%Nm/rad
wn=sqrt(k/m);
s_=0.25;
%wd=sqrt(1-s_^2)*wn;
c=2*s_*wn*m;
%系统方程:
H=[0,1,0,0];
%Covariances of w:
Q=diag([0,Ts*0.1*Ts',0,0]);
%Covariances of v:
R=[0.1];
%初始估计值:
x=[0;0;k/2;c];
%初始估计误差协方差:
P=diag([0,Ts*0.1*Ts',40000,0.000000]);
u=x(2:4);%开始时刻不进行卡尔曼滤波,直接将初始值输出,2020
end
%%卡尔曼滤波::
if t>0
%先验,Time update:
%根据系统状态方程计算下一状态预测值:
%x=A*x+B*u1;
Jf=eye(4);
Jf(1,2)=Ts;
Jf(2,1)=-x(3)/m*Ts;
Jf(2,2)=1-x(4)/m*Ts;
%*********very important**********
%因为激励为飞轮的转角与转速而不是飞轮对离合器的激励转矩,改动如下:
%Jf(2,3)=-x(1)/m*Ts;
Jf(2,3)=(s_flywheel-x(1))/m*Ts;
%Jf(2,4)=-x(2)/m*Ts;
Jf(2,4)=(v_flywheel-x(2))/m*Ts;
%**********************************
s=x(1);
v=x(2);
k=x(3);
c=x(4);
F_input=c*v_flywheel+k*s_flywheel;
s=s+Ts*v;
v=-k/m*Ts*s+(1-c/m*Ts)*v+Ts/m*F_input;
x(1)=s;
x(2)=v;
x(3)=k;
x(4)=c;
%************强跟踪滤波器STF核心代码《::*********************************
%**************************************
% 可调节参数汇总
%**************************************
rho=0.95; % 残差方差阵的遗忘因子初始值 (0.95<rho<0.995)
beta=1; % 量测噪声的衰减因子选定值 (beta>=1)
%beta_upmax=8; % 衰减因子选定极限
%***************************************
F=Jf;
r0=z-H*x;
if t<2*Ts
V0=r0*r0';
else
V0=(rho*V0+r0*r0')/(1+rho);
end
N=V0-beta*R-H*Q*H';
M=H*F*P*F'*H';
eta=trace(N)/trace(M);
if eta>1
lamda=eta;
else
lamda=1;
end
%预测误差协方差:
P=lamda*Jf*P*Jf'+Q;
%!!!渐消因子lamda=1即为普通卡尔曼滤波!!!
%************STF核心代码》************************************************
%后验,Measurement update:
%根据估计误差协方差和测量噪声协方差计算卡尔曼增益:
Kk=P*H'/(H*P*H'+R);
%计算最优估计值:
%x=A*x+Mn*(yv-C*A*x);
x=x+Kk*(z-H*x);
%计算最优估计值和真实值之间的误差协方差矩阵,为下次递推做准备:
P=(eye(4)-Kk*H)*P;
%ye=C*x; %Filtered value
u=x(2:4); %Filtered signal
errcov=H*P*H'; %Covariance of estimation error
end
%}