智能车制作

标题: 授之以渔： 卡尔曼滤波器 ...大泻蜜 ... [打印本页]
作者: Oner    时间: 2011-12-14 14:37
标题: 授之以渔： 卡尔曼滤波器 ...大泻蜜 ...
本帖最后由 Oner 于 2011-12-14 14:41 编辑

首先声明：原创作者是highgear，独乐乐不如众乐乐，于是乎，看到好东西当然要和大家分享啦。
希望对各位电磁组的大大们有些许帮助咯。
作者: Oner    时间: 2011-12-14 14:37
本帖最后由 Oner 于 2011-12-14 14:39 编辑

（二）
一片绿油油的草地上有一条曲折的小径，通向一棵大树。一个要求被提出：从起点沿着小径走到树下。

“很简单。” A说，于是他丝毫不差地沿着小径走到了树下。

现在，难度被增加了：蒙上眼。

“也不难，我当过特种兵。” B说，于是他歪歪扭扭地走到了树 ………. 旁。“唉，好久不练，生疏了。”

“看我的，我有 DIY 的 GPS！” C说，于是他像个醉汉似地走到了树………. 旁。“唉，这个 GPS 软件没做好，漂移太大。”

“我来试试。” 旁边一人拿过 GPS,  蒙上眼，居然沿着小径走到了树下。

“这么厉害！你是什么人?”

“卡尔曼 ! ”

“卡尔曼？！你是卡尔曼？”众人大吃一惊。

“我是说这个 GPS 卡而慢。”
作者: Oner    时间: 2011-12-14 14:38
（三）
首先回顾一下传统数字滤波器。

对于一个线性时不变系统，施加一个输入 u(t) ，我们可以得到一个输出 y(t) . 如果输入是一个冲击，则输出y(t) 被称作冲击响应，用 h(t) 来表示，是系统的内核。对于任意 u(t), 输出 y(t) 可以通过 u(t) 与冲击响应 h(t) 的卷积得到，这是 FIR 滤波器的基本原理。我们还可以通过系统微分方程转换为差分方程，或是通过 laplace 传递函数转换到差分方程，最后得到一个递推公式，这种形式的滤波器就是IIR 滤波器。
以前讲过，一个系统可以用时域的微分方程来建立，然后可以用 laplace 的传递函数来处理，把解微分方程变为多项式乘法，可以简单的求解。还有另外一种处理形式就是状态空间，以矩阵形式来处理微分方程或微分方程组，利用矩阵变换求解，类同齐次方程组的矩阵形式。例如微分方程:

              y’’ +3y’ + 2y = u

让      X1 = y,       X2 = y’ = X1’,   则上式变为：
             X2’ = -3 X2  – 2 X1  - u
             X1’ =      X2
矩阵形式为：


通用形式为：
           X’ = A*X + B*u
           Y = C*X.

可以看到，可以很轻易的微分方程或微分方程组转换到状态空间形式，而状态空间与laplace 传递函数之间可以相互转换，事实上
矩阵A 的特征值就是s传递函数的极点。系统的传递函数（阵）可以通过矩阵变换得到：

               Y(s) = C * (s * I - A) -1 * B

同理，连续域的微分方程对应了离散域的差分方程，s 对应了z,   离散域状态空间相应的变为：
              X(k) = A*X(k-1) + B(u-1)
              Y(k) = C*X(k)
作者: Oner    时间: 2011-12-14 14:38
（三）
我们现在来看看蒙眼走小径的走法问题。


假设A 走过的路径是真真正的路径，为Za; B是用自己的大脑作为预测估计器，走出了一个预测路径，
为 Zb; C 用测量器，走出了一条测量路径，为Zc。用图片来说明：




“系统真实输出”是 A 走过的路径： Za = C * X；
“测量输出”是Zb.   Zb = Za + V，
这里 V 是噪声，即GPS 的漂移；
“预测估计输出”是 Zc = C * X^,
X^是预测的状态。
T 是采样延时。

现在，蒙上眼的情况下有两种选择，GPS 或大脑预测估计器。如果GPS很准而预测不准，那么可以选择GPS；
如果预测准确而GPS不准，那么选择预测估计器， 等等， 没有反馈的预测估计器会因为累积误差而导致越来越不准。如果两个都不准，该如何取舍？如何把两者结合在一起呢？

我们可以设置一个信心指数 K，K 在 0 与1之间，来说明对测量值还是预测值的信任程度：

        Z = K * Zb + (1 – K) * Zc  = Zc + K*(Zb – Zc)           (1)

可以看出， 当 K = 1 和 0 时，分别选择了GPS 或预测估计器. 现在，可以把误差 Zb -Zc 作为反馈误差，来修正 预测估计器的结果。
新的系统结构图如下：




这个框图，就是卡尔曼滤波器的基本构造。学过现代控制理论的同学都这个图应该很熟悉，与状态变量估计控制的图形差不多，只是其中的 K = 1 而且没有噪声项和系统反馈而已。
而我们下面的任务，就是如何确定这个K 值。
作者: perfect_co    时间: 2011-12-14 14:39
看不懂
作者: Oner    时间: 2011-12-14 14:39
本帖最后由 Oner 于 2011-12-14 14:42 编辑

关于卡尔曼滤波器的推导过程，枯燥晦涩，我就略过，直接关注结果。




（四）
计算过程：

卡尔曼滤波是一种递归的估计，即只要获知上一时刻状态的估计值以及当前状态的观测值就可以计算出当前状态的估计值，因此不需要记录观测或者估计的历史信息。


卡尔曼滤波器的递归过程：

1)   估计时刻k 的状态：
          X(k) = A*X(k-1) + B*u(k)
     这里,   u(k) 是系统输入

2)   计算误差相关矩阵P, 度量估计值的精确程度：
          P(k) = A*P(k-1)*A’+ Q
      这里,   Q = E{ Wj^2 } 是系统噪声的协方差阵，即系统框图中的Wj的协方差阵, Q 应该是不断变化的，为了简化，当作一个常数矩阵。

3)   计算卡尔曼增益, 以下略去 (k),  即 P = P(k), X = X(k)：
             K = P *C’ * (C * P * C’ + R) -1
       这里 R = E{ Vj^2 }, 是测量噪声的协方差(阵),  即系统框图中的 Vj 的协方差, 为了简化，也当作一个常数矩阵。由于我们的系统一般是单输入单输出，所以 R是一个 1x1的矩阵，即一个常数，上面的公式可以简化为:
             K = P *C’ / (C * P * C’ + R)

4)     状态变量反馈的误差量：
             e = Z(k) – C*X(k)
         这里的 Z(k) 是带噪声的测量量

5)    更新误差相关矩阵P
          P = P – K * C * P

6)   更新状态变量：
            X =X + K*e = X + K* (Z(k) – C*X(k))

7)   最后的输出：
            Y = C*X


现在的问题就是如何实现卡尔曼滤波， A, B, C, Q, R 这些矩阵或量如何确定？
作者: Oner    时间: 2011-12-14 14:40
本帖最后由 Oner 于 2011-12-14 14:42 编辑

（五）
仿真实例

下面用仿真实例来观察卡尔曼滤波器的效果。假设我们的系统是一个加热系统，热时间常数为 60 秒，100度时达到热平衡。忽略系统的延迟，那么当系统加电后，温度由 0 开始上升。这个上升过程大家应该很熟悉，这是一个指数函数：

        y(t) = 100 * (1 – e-t/60)

其 laplace 传递函数为:
      
        y(s) = 100 / (60 * s - 1)

我们人为的加入了随机噪声来模拟测量噪声





我们假定并不知道系统的传递函数，现在只是简单的, 随便地构造了一个预测系统。
A = [1, 0; 0, 1]
B = [1; 0]
C = [1, 0]

这是一个二阶系统，其输出是一条直线，与实际的系统相去甚远：




测量噪声的协方差 R = 40，此为猜测值; 系统噪声 Q = 2，也是猜测值，预测模型越不准，Q 值应越大。
卡尔曼滤波器的结果，红色为滤波器输出：



可以看到，尽管我们使用了一个粗劣的预测估计器，Kalman 滤波器还是相当的皮实，基本上消除了噪声.
作者: Oner    时间: 2011-12-14 14:40
精确模型的建立
要建立一个精确的预测估计模型，我们还是要利用方差。如果一个估计的曲线与实际曲线完全重合时，他们的方差为 0. 方差越小，拟合度越高， 最小二乘法的原理便是如此。具体推导过程还是省略，直接给出 matlab 的拟合程序，这是一个非常非常有用的程序。
如果数学模型很精确， 能不能直接数学模型的输出作为滤波器的结果呢？不能，因为没有反馈，数学模型的输出会因为没有反馈的校正造成误差不断累积，失之毫厘，谬之千里。

下面是用最小二乘法获得系统的模型并做为预测估计器，设定为 3 阶系统， 得出的数学模型相当准确，所以Q值可以取一个小值，
这里 Q= 0.02, 现在看看卡尔曼滤波器的结果：




效果非常好，卡尔曼滤波器的输出与实际系统的输出 (即无噪声的系统输出) 几乎重合，这是精确的预测估计模型带来的好处。

现在比较两个例子中 卡尔曼增益的不同



最小二乘法获得系统的模型中的增益迅速地由大变小，最后小于不准确模型。K 值较小，意味着误差反馈量较小， 使得预测输出更偏重信任预测模型的结果， 通过上图可以看到 kalman 算法的自适应性。至于误差相关阵 P 的值，同样，最小二乘法获得系统的模型中的 C*P*C' 的值较小， 这里就不给出图形了。

结论：
从上面的例子可以看出，卡尔曼滤波器对于预测系统的要求并不高，所以，那个”卡尔曼”可以蒙上眼，拿着一个劣质的 GPS 可以走到树下。若系统的预测模型很准确，卡尔曼滤波器会有一个相当良好的效果。
作者: Oner    时间: 2011-12-14 14:42
占个座。以后可能有用。
作者: xjtuzhanghongji    时间: 2011-12-15 09:32
本帖最后由 xjtuzhanghongji 于 2011-12-15 09:34 编辑

写的确实很不错，顺便说下，前半部分公式里面的’应该导数的意思，包括两个’’的二阶导数意思，而后半部分的公式里的’是向量转置的意思，大家不要搞混了
而-1应该是上标，是逆矩阵的意思吧，如果考虑到某些矩阵实际是1维的（标量），所以-1变成了除号
作者: chen19910528    时间: 2011-12-15 15:12
汗！我看了几个小时才把卡尔曼滤波器弄懂！！
作者: whut_RY    时间: 2011-12-15 16:10
MARK!
作者: skul_malke    时间: 2011-12-15 17:16
学习了
作者: zl6977    时间: 2011-12-16 15:35
学习了
作者: miss126    时间: 2011-12-16 17:30
回复 9# Oner


    没有图片啊 不错写的确实不错！！ 我大概理解了
作者: 朱杰亮    时间: 2011-12-16 18:40
不懂 回去好好看看
作者: wjyeasy    时间: 2011-12-16 20:44
这内涵贴看的费劲啊
作者: liughsy    时间: 2011-12-16 23:10
此位有人占座，谢谢合作！
作者: liughsy    时间: 2011-12-16 23:10
此位有人占座，谢谢合作！
作者: zhu_xuekui    时间: 2011-12-16 23:11
求图片看看实际效果。。。楼主强悍。。。佩服  佩服
作者: 594961689    时间: 2011-12-17 02:20
顶！！！
作者: 7998    时间: 2011-12-17 12:00
这个，，，看起来挺费劲的
作者: zouyf12    时间: 2011-12-18 18:18
卡尔曼原来是需要建模的啊，  这个很麻烦啊
作者: hanahimiEX    时间: 2011-12-18 19:20
多谢啦
作者: hanahimiEX    时间: 2011-12-18 19:20
谢谢啦
作者: 我是传奇    时间: 2011-12-18 19:25
matlab是个好东西，建议兄弟们好好学习下
作者: 晨景的风    时间: 2011-12-18 20:17
好吧，看不懂
作者: aaaalook    时间: 2011-12-19 21:35
mark一下，哈哈，很高深啊。
http://bbs.21ic.com/icview-292853-1-1.html
作者: aaaalook    时间: 2011-12-19 21:37
有木有？
第一个例子中的代码
  1. Ts = 1;                    采样时间
  2. s = tf('s');
  3. sysc = 100/(60*s +1);      真实系统的传递函数
  4. sysd = c2d(sysc, Ts);

  6. t = 0:TsTs*300);
  7. u = ones(1, length(t));     系统输入
  8. ys = lsim(sysd, u, t, 0);   系统输出
  9. ys = ys + 10*rand(size(ys));  测量量 = 系统输出 + 噪声

  11. A = [1, 0; 0, 1];
  12. B = [1; 0];
  13. C = [1, 0];
  14. N = length(A);          系统维数
  15. Len = length(u);

  17. yout = zeros(Len, 1);   滤波器输出
  18. Xh = zeros(N, 1);       状态变量
  19. P = eye(N);              
  20. Q = 2 * eye(N);         系统噪声
  21. R = 50;                 测量噪声
  22. for i = 1 : Len
  23.     Xh = A*Xh + B*u(i);
  24.     P = A*P*A' + Q;
  25.     K = P*C'* inv(C*P*C' + R);
  26.     Xh = Xh + K*(ys(i) - C*Xh);
  27.     P = P - K*C*P;
  28.     yout(i) = C*Xh;
  29. end
  30. Ts = 1;                    采样时间
  31. s = tf('s');
  32. sysc = 100/(60*s +1);      真实系统的传递函数
  33. sysd = c2d(sysc, Ts);

  35. t = 0:Ts:(Ts*300);
  36. u = ones(1, length(t));     系统输入
  37. ys = lsim(sysd, u, t, 0);   系统输出
  38. ys = ys + 10*rand(size(ys));  测量量 = 系统输出 + 噪声

  40. A = [1, 0; 0, 1];
  41. B = [1; 0];
  42. C = [1, 0];
  43. N = length(A);          系统维数
  44. Len = length(u);

  46. yout = zeros(Len, 1);   滤波器输出
  47. Xh = zeros(N, 1);       状态变量
  48. P = eye(N);              
  49. Q = 2 * eye(N);         系统噪声
  50. R = 50;                 测量噪声
  51. for i = 1 : Len
  52.     Xh = A*Xh + B*u(i);
  53.     P = A*P*A' + Q;
  54.     K = P*C'* inv(C*P*C' + R);
  55.     Xh = Xh + K*(ys(i) - C*Xh);
  56.     P = P - K*C*P;
  57.     yout(i) = C*Xh;
  58. end

  60. plot(t, ys, t, yout);
复制代码
作者: aaaalook    时间: 2011-12-19 21:38
使用最小二乘法例子的代码：


  3.     Ts = 1;               采样时间
  4.     s = tf('s');
  5.     sysc = 100/(60*s +1);   真实系统的传递函数
  6.     sysd = c2d(sysc, Ts);

  8.     t = 0:Ts:(Ts*300);
  9.     u = ones(1, length(t));   系统输入
  10.     ys = lsim(sysd, u, t, 0); 系统输出
  11.     ys = ys + 10*rand(size(ys)); 测量量 = 系统输出 + 噪声

  13.     [numd, dend] = LeastSquare(u, ys, 3);   最小二乘法获取预测系统模型
  14.     [A, B, C, D] = tf2ss(numd, dend);       变换到状态空间形式
  15.     Len = length(u);
  16.     N = length(A);         系统维数

  18.     yout = zeros(Len, 1);  滤波器输出
  19.     Xh = zeros(N, 1);      状态变量
  20.     P = eye(N);
  21.     Q = 0.02 * eye(N);     系统噪声
  22.     R = 50;                测量噪声
  23.     for i = 1 : Len
  24.         Xh = A*Xh + B*u(i);
  25.         P = A*P*A' + Q;

  27.         K = P*C'* inv(C*P*C' + R);
  28.         Xh = Xh + K*(ys(i) - C*Xh);
  29.         P = P - K*C*P;
  30.         yout(i) = C*Xh;
  31.     end

  33.     plot(t, ys, t, yout);
复制代码
作者: aaaalook    时间: 2011-12-19 21:40
最小二乘法获取系统传递函数：
  1. function [numd, dend] = LeastSquare(x, y, N)
  2. count = length(y);
  3. M = count - 1;
  4. ai = zeros(N*2, M);
  5. for i=1:N
  6.     ai(i, i:M) = x(1:(count-i));
  7.     ai(i+N, i:M) = -y(1:(count-i));
  8. end
  9. bi = y(2:(count));
  10. xd = (inv(ai*ai')*ai*bi)';
  11. numd = [0 xd(1:N)];
  12. dend = [1 xd(N+1: N+N)];
复制代码
输入参数：
x  ---  系统输入阵列
y  ---  系统输出阵列
N  ---  系统传递函数的阶数

函数输出：
系统 Z 传递函数分子与分母的系数，即： h(z) = numd(z) / dend(z)

这个函数实在是居家旅行 ..., 啊，是建模仿真必备良器。
作者: 篪骋    时间: 2011-12-20 00:36
作者: 无与伦比119    时间: 2011-12-20 15:55
以后肯定有用的，mark了
作者: 虫总    时间: 2011-12-20 17:15
无图无真相
作者: wenhaoyuan    时间: 2011-12-22 23:07
强贴留名
作者: 邛于    时间: 2012-2-14 15:08
学习中。。。。。。
作者: 断翅at雄鹰    时间: 2012-2-14 21:05
:victory::Q:Q
作者: 津生有你    时间: 2012-2-15 11:14
作者: fanyu19900706    时间: 2012-2-15 13:21
楼主，大神啊
作者: 我の小车    时间: 2012-2-15 17:50
顶下，谢了啊
作者: yu2010550221    时间: 2012-2-25 19:18
作者: 杜道轶    时间: 2012-4-10 11:23
作者: zhuwenwujy    时间: 2012-4-10 12:27
高端……MARK
作者: sweetgum    时间: 2012-4-14 12:44
Mark
作者: ou421293238    时间: 2012-4-21 09:35
有没有有关最小二乘法在转向舵机中的应用呢？
作者: 白蓝鸽    时间: 2012-4-21 20:33
感谢楼主这么用心， 虽然没看懂（因为基础太差，正在学） 还是要吼吼
作者: sunnyfei    时间: 2012-6-26 09:58
学习
作者: gggfff    时间: 2012-6-26 11:21
略懂略懂



欢迎光临 智能车制作 (http://111.231.132.190/) Powered by Discuz! X3.2