求助：有关扩展卡尔曼滤波辨识的收敛速度问题

woshi217

是这样的，我用S函数编写了一个扩展卡尔曼滤波辨识电机转动惯量的程序，辨识结果收敛时间大概在7-8s的样子，但是看到国内外很多论文上大概的收敛时间在1s以内，个人感觉问题所在已在参考程序中注明，但是没有想到解决方法，所以希望知道的朋友能够帮忙解决一下，先谢谢了！




参考程序：function [ sys, x0, str, ts ] = EKFLITER(t, x, u, flag)
%This S - Function is to indentify the inertia online using extended kalman
%filter method.
%x(1) = w(estimate speed)
%x(2) = TL(load torque)
%x(3) = 1/J(J is inertia of motor)
%u(1) = Te
%u(2) = w(the real speed)
switch flag,
    case 0,
        [ sys, x0, str, ts ] = mdlInitializeSizes;
    case 2,
        sys = mdlUpdate(t, x, u);
    case 3,
        sys = mdlOutputs(t, x, u);
    case {1, 4, 9},
        sys = [ ];
    otherwise
        error([ 'Unhandled flag = ', num2str(flag) ]);
end
function [ sys, x0, str, ts ] = mdlInitializeSizes
sizes = simsizes;
sizes.NumContStates = 0;
sizes.NumDiscStates = 3;
sizes.NumOutputs = 2;
sizes.NumInputs = 2;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;
sys = simsizes(sizes);
x0 = [ 0 0 1500 ];
str = [ ];
ts = 1e-4;
function sys = mdlUpdate(t, x, u)
ts = 1e-4;
A = [ 1 - x(3)*ts 0; 0 1 0; 0 0 1 ];
B = [ x(3)*ts 0 0 ]';
P = diag([ 1 1 10 ]);                                    %*
C = [ 1 0 0 ];
H = [ 1 0 0 ];
Q = diag([ 1e-5 1e-4 5e-6 ]);
R = 5e-3;
y = u(2);
Y = C*x;
G = [ 1 - x(3)*ts u(1) - x(2)*ts; 0 1 0; 0 0 1 ];
%The first step
x = A*x + B*u(1);

P = G*P*G' + Q;    %*

% P = G*P*G' + Q;%这里应该是用P = G*Pe*G’ + Q但是由于Pe在之前没有定义所以不能用, 现在的问题在于, 本来是一个迭代循环, 可是现在P每次更新都被重置为diag[ 1 1 10 ], 所以程序的收敛速度非常慢。看看哪位朋友能帮忙解决一下, 就是起始的P值为diag[ 1 1 10 ], 之后P = G*Pe*G’ + Q, 多谢大家了, 已经弄了好几天了, 没有找到方法。
% %The second step

K = P*H'/(H*P*H' + R);
Xe = x + K*( y - Y );
Pe = P - K*H*P;
%Result
sys(1) = Xe(1);
sys(2) = Xe(2);
sys(3) = Xe(3);
function sys = mdlOutputs(t, Xe, u)
%Output sys(1) is load torque
%sys(2) is inertia J
sys(1) = Xe(2);
sys(2) = 1/Xe(3);

Rainyboy

回复 1 # woshi217 的帖子

Pe是P迭代到下一步之后的值么？
我是说，为什么楼主说“Pe不能定义”呢？

不了解simu，还请楼主解释一下，共同学习

woshi217

回复 2 # Rainyboy 的帖子

P是先验估计误差的协方差矩阵，Pe是后验估计误差的协方差矩阵，两者是不同的；
第一步运算（时间更新）中的Pe是预估计值，也就是给定值，而第二部运算（状态更新）中的Pe是由第一步中的P值通过公式运算来得到的，而且整个迭代过程中只需要第一次运算给定Pe值就行了，之后的运算都是由公式P - K*H*P得到的，但是整个function sys = mdlUpdate(t, x, u)本身就是一个循环过程，所以不知道Pe的第一次初值应该怎么添加进来？（以上观点仅代表个人想法，如有错误请指出，谢谢！）

Rainyboy

回复 3 # woshi217 的帖子

是不是说，是这样的问题，下一步计算P的时候需要上一步计算P的结果，但是现在不知道怎么实现，就导致每次计算的时候P都是重新初始化了的？
本来应该是：
Pint->p1->P2->p3...
现在却是：
Pint->p1
Pint->p2
是么？

woshi217

回复 4 # Rainyboy 的帖子

是的，可以这么理解，您有什么好的处理方法吗？

Rainyboy

回复 5 # woshi217 的帖子

你所用的环境中能使用全局变量么？
如果可以的话，声明一个全局变量，储存上一步计算之后的结果不就好了？
如果没有，将上一步的计算结果返回，例如将函数的接口成这样：
[sys P_new] = mdlUpdate(t, x, u, P_last)

那么，在初始化的时候：
    P_last = 初始值；
    [sys,P_new] = mdlUpdate(t, x, u, P_last)；
    P_last = P_new;
在后续的调用中。不执行初始化的那句话，加一个条件判断就可以吧。只执行：
    [sys,P_new] = mdlUpdate(t, x, u, P_last);
    P_last = P_new;

是不是就可以了？