Fundmentals of Kalman Filtering 翻译1

Chapter 1

介绍

      本章我们将致力于讨论本书中涉及的所有数学理论。我们首先介绍如何将两种不同的数字积分技术应用于线性和非线性微分方程。当我们必须积分用来评估了卡尔曼滤波的微分方程时，数字积分技术非常必要。另外，从非线性微分方程中得到状态量有时需要用到数字积分技术。接下来，介绍状态空间符号用来表示不同的微分方程。在卡尔曼滤波之前，现实世界的问题必须转化成状态空间符号。最后，给出计算基础矩阵的方法。基础矩阵可以将卡尔曼滤波转化成数字形式。

微分方程的数值积分

      本书中我们将会模拟线性和非线性常微分方程。由于这些方程没有封闭解，所以有必要求助于数值积分技术去解或者模拟这些方程。许多数值积分技术就是为了解微分方程而存在，让我们从最简单的欧拉积分开始。

      欧拉积分是一种可以应用于本书中绝大多数例子的数值积分技术。考虑一个一阶微分方程：

 。我们从函数的导数定义中可以知道：

 ，h 为 很小的数。

重写上式得到著名的欧拉积分方程： 

。该方程描述了在下一积分间隔，状态量 

 的值，与之前时刻的状态量加之前时刻导数值的比例 有关。每次循环，时间加上积分间隔 h 。为了开始积分过程，我们需要知道 

 的初始值。初始值从微分方程得到。

      我们知道，对一个函数积分相当于，在坐标系中求它的曲线下方的面积。欧拉积分是计算

 所代表的曲线下方面积的最简单的近似方法。用高为

，宽为 h 的矩形来表示区域的面积，如图1.1所示。

Fig 1.1

      当然用矩形近似曲线面积会带来误差，但是误差可以通过减小 h 来减小。减小h 相当于让积分步长更小。

      积分步长必须相当小，来保证计算结果满足要求的计算精度。在工程实践中最简单的设置积分步长的方法是进行实验测试。根据经验法则，初始步长一般比考虑中的系统中最小的时间常数小几倍。然后将步长减小一半，看结果是否有明显变化，如果结果近似相等，那么采用更大的步长来避免过度的计算时间。如果结果提升明显，就将积分步长在此基础上再减小一半，重复整个过程。

      为了直观感受欧拉积分如何解微分方程，让我们来解决一个已经知道答案的问题。如果关于 x 的微分方程的解由正弦曲线给出： 

  ，则微分方程为： 

 ，在次微分得到 ：

 ，对比可得线性二阶微分方程：

，容易求得初始状态：

 。有了这系已知条件，让我们看看数字积分是否能得到正确的结果。

      为了检验 

 的封闭解的理论，给出了一个关于数字积分的模拟程序（MATLAB代码）。一个使用欧拉数字积分技术的二阶微分方程的模拟程序由Listing 1.1 给出。首先，表示两阶微分方程。我们先求解 

 的一阶导数，然后求解 

 。然后时间加上积分间隔 

 ，我们重复程直到时间大于10s。我们每隔0.1s 同时输出理论和模拟值（欧拉积分得到）。

%Listing 1.1 
%simulation of second-order differential eqution using euler numerical
%intergration
clear all
close all
count = 0;
w=2;
t=0;
s=0;x=0;
xd=w;
h=0.01;
while t<=10
   s=s+h;
   xdd=-w*w*x;
   xd=xd+h*xdd;
   x=x+h*xd;
   t=t+h;
   if s>=.09999
      s=0;
      xtheory = sin(w*t);
      count=count+1;
      ArrayT(count)=t;
      ArrayX(count)=x;
      ArrayXtheory(count)=xtheory;  
   end  
end
figure
plot(ArrayT,ArrayX,ArrayT,ArrayXtheory),grid
xlabel('time(sec)')
ylabel('x&xtheory')
axis([0,10,-2,2])
clc
output = [ArrayT',ArrayX',ArrayT',ArrayXtheory'];
save datfilL1.txt output -ascii 
disp(count)
disp 'simulation finished'
           

      上面代码运行结果如下：

       可以看出来，积分步长 

=0.01 时，欧拉数字积分准确的积分了二阶微分方程，这是由于精确解和数值解是一致的。像之前提到的，由于我们经常没有精确解，所以通常通过实验确定积分步长。

      在给定积分间隔时，对比其他数值积分方法，有些时候使用欧拉积分的代价是精度较低。当使用欧拉积分时，谨慎并使用小的积分步长是提高精度的答案。然而，有些时候不论积分步长多小，欧拉积分完全不能使用。由于这个原因我们考虑使用精度更高的数值积分技术 二阶龙格库塔 法。

      两阶龙格库塔数值积分过程可以很简单的进行状态化。给出一阶微分方程：

 ，t 是时间，我们试图寻找关于时间的函数x 的递归关系。根据二阶龙格库塔数值积分技术，x 在有、下个积分间隔h 的值由下式给出：

下标k-1 代表上一时刻，k 代表当前时刻。通过上式我们可以知道，当前时刻x 值等于上一时刻x 值加上一定比例的  t 时刻和 t+h 时刻 的导数值。为了开始积分，我们必须知道 x 的初始值。就像之前提到的，初始值取决于微分方程的初始状态。

      如果观察图1.3， 我们可以发现二阶龙格库塔数值积分其实也是一种寻找曲线下的面积。然而，我们使用小梯形而不是矩形来近似。因为我们可以将梯形的面积看做时间为底，与两个高之和的面积的一半，我们可以看出来二阶龙格库塔法的积分形式实际上是梯形的面积公式。图1.3 也表明了用梯形近似曲线之下的面积胜过用矩形近似。因此，对于给定的积分步长，二阶龙格库塔法比欧拉积分精度更高。

      一个使用为而阶龙格库塔数值积分的仿真列表在Listing 1.2 给出。从仿真程序中可以看出微分方程和导数信息出现了两次并用粗体标出。在一个积分间隔我们计算了微分方程两次：一次计算t 时刻的导数，一次计算t+h 时刻的导数。可以看出，程序中每隔0.1秒输出x的微分方程的数直结果利用变量XTHEORY。在这个例子中，发微分方程中的自然频率w值为 2rad/s(W=2)。 

%Listing 1.2 
%simulation of second-order differential eqution using runge-kutta
%intergration
clear all
close all
count = 0;
w=2;
t=0;
s=0;x=0;
xd=w;
h=0.01;
while t<=10
   s=s+h;
   xold=x;
   xdold=xd;
   xdd=-w*w*x;
   x=x+h*xd;
   xd=xd+h*xdd;
   t=t+h;
   xdd=-w*w*x;
   x=0.5*(xold+x+h*xd);
   xd=0.5*(xdold+xd+h*xdd);
   if s>=.09999
      s=0;
      xtheory = sin(w*t);
      count=count+1;
      ArrayT(count)=t;
      ArrayX(count)=x;
      ArrayXtheory(count)=xtheory;  
   end  
end
figure
plot(ArrayT,ArrayX,ArrayT,ArrayXtheory),grid
xlabel('time(sec)')
ylabel('x&xtheory')
axis([0,10,-2,2])
clc
output = [ArrayT',ArrayX',ArrayT',ArrayXtheory'];
save ('datfileC1L2.txt' ,'output', '-ascii')
disp(count)
disp 'simulation finished'
           

  运行Listing1.2 ，结果如图1.4 所示，从结果可以看出：给定积分步长（H=0.01），由于真值和积分结果一致，两阶龙格库塔数值积分方法准确的积分了二阶微分方程。

      我们将在正本书中使用两阶龙格库塔法，因为它容易理解、编程方便，最重要的是对本书中所有的例子都能产生准确的结果。所以，本书中后续涉及数值积分的程序将会与Listing 1.2 有相同的结构。

状态空间符号

基础矩阵

总结

相关文献

待续。。。。。