《确定高精度参数问题论文设计.doc》由会员分享，可在线阅读，更多相关《确定高精度参数问题论文设计.doc（20页珍藏版）》请在三一文库上搜索。
1、确定高精度参数问题摘要本文主要研究了在要求精确制导等有关高科技的实际研究背景中，如何建立数学模型以及如何高精度地估计数学参数问题。问题一，在观测数据无误差的情况下，将参数估计问题转化为最小二乘求解微分方程组，得到两生物数目的隐式方程式。得到待定参数的值为，,，。问题二，在观测数据无误差的情况下,4个参数的值均未知，利用问题一的解决方法，将4个参数的确定问题转化为1个比例系数C的求解问题。然后用龙格库塔法进行数值积分确定出比例系数应满足的条件，并利用搜索寻优的方法进行求解，得到待定参数的值为，。问题三，针对观测时间准确而观测资料有误的特点，建立基于数据去噪，将参数识辨问题转换为最小方差意义下的参
2、数估计问题。通过对该模型的分组平滑改进，高精度的估计出参数的数值。进而采用固定点最优平滑方法，从而得到参数的数值。利用DATA2得到最终的参数估计值为：，。问题四，问题中观测时间和观测资料均存在误差的限制条件，建立包含观测时间状态变量的推广型卡尔曼滤波模型，用第三问题中改进的模型，计算出参数的数值。采用固定点最优平滑方法，确定出初始值，从而得出改组观测数据下，最小方差意义的系统模型参数的数值。最终参数估计结果为：，。关键字：参数估计、最小二乘、卡尔曼滤波模型、最优平滑1 问题重述考虑航天器在仅受到地球万有引力、航天器自身发动机作用力的作用下作平面运动，将地球和航天器视为质点，由理论力学可知，一
3、个刚体在空间的运动可以看作质心的移动,因此可以应用质心运动定理来研究刚体质心的移动规律。以地球中心为原点，建立直角坐标系，航天器绕地球飞行，可以出现在该直角坐标系中四个象限的任意一个之内。平面直角坐标系如图1。符号说明如下:航天器在x方向的速度航天器在y方向的速度万有引力，航天器发动机作用力，为控制变量万有引力与x轴正方向的夹角航天器发动机作用力与x轴正方向的夹角初始时刻，航天器初始位置航天器x方向初速度航天器y方向初速度航天器受的万有引力方向指向地球中心（原点），航天器受推力的方向与x轴正方向成角。将和投影到该直角坐标系上，见图1图1 航天器受力分解图其中， 初始条件为，都是关于时间的位置函
4、数，航天器在x方向的分速度即对时间t求导，航天器在y方向的分速度即对时间t求导，航天器在x方向的加速度即对时间t求二阶导，航天器在y方向的加速度即对时间t求二阶导。由此建立的航天器模型如下：显然这样的数学模型在精度上是远远不能满足实际需要的，在其他要求精确制导等有关高科技的实际问题中，我们都面临着类似的问题：我们必须建立高精度的数学模型，必须高精度地估计模型中的大批参数，因为只有这样的数学模型才能解决实际问题，而不会出现差之毫厘，结果却失之千里的情况。这时所建立数学模型的精度就成了数学模型的生命线。例如上述问题中的航天器还要受到地球质量分布不均匀所引起的摄动力，大气阻力，日、月及其它星球的摄动
5、引力的影响，以及航天器发动机为调整航天器自身姿态运作时作用力的影响。这样不但数学模型十分复杂，而且在这些数学模型中还要涉及到许多重要的参数，如地球的引力场模型就有许多待定参数。不仅如此，在对航天器进行测量时，还涉及到观测站的地理位置以及设备的系统误差等参数。为此人们要设法利用长期积累的丰富的观测资料，高精度确定这些重要的参数。由于航天器的问题太复杂，本题仅考虑了一个较简单的确定高精度参数问题。假设有一个生态系统，其中含有两种生物，即： A生物和B生物，其中A生物是捕食者，B生物是被捕食者。假设时刻捕食者A的数目为，被捕食者B数目为，它们之间满足以下变化规律：初始条件为： ,其中为模型的待定参数
6、。通过对此生态系统的观测，得到了相关的观测数据。现要求解决如下问题：1) 在观测数据无误差的情况下，若已知，求其它5个参数。有关数据见数据文件：DATA1.TXT2) 在观测数据无误差的情况下，若也未知，问至少需要多少组观测数据，才能确定参数。有关数据见数据文件：DATA1.TXT3) 在观测资料有误差（时间变量不含有误差）的情况下，请分别利用观测数据DATA2.TXT和DATA3.TXT，确定参数在某种意义下的最优解，并与仿真结果比较，进而改进所建的数学模型。4) 假设连观测资料的时间变量也含有误差，利用数据DATA4.TXT，建立数学模型，确定参数在某种意义下的最优解。2 模型假定虽然该问
7、题已经给定食饵捕食者的具体模型，但是为了问题分析与求解的方便我们做出如下假定：1、假定DATA1所给数据无误差；2、假定DATA2、DATA3、DATA4所给数据仅含有随机误差而不含系统误差，且其统计特性与高斯白噪声特性相同；3、非线性微分方程的理论解一定存在，而且这个理论解与实际解之间能够用一个线性微分方程表示，此时我们可以说，理论解能够“充分地”对系统的实际特性给予分析。3 符号定义：捕食者、食饵分别独立生存时的自然增长率。：食饵对捕食者的供养能力。：捕食者掠取食饵的能力。：捕食者、食饵在观测初始时刻的数目。：分别表示捕食者、食饵的数目。：在本模型中表示待定参数，状态变量和设计各部分各因素
8、的索引下标。：该食饵捕食者生物系统的状态变量与输出变量。：系统噪声与观测噪声。：系统状态变量对时间导数的连续型非线性函数。X、Y：该生物系统的状态向量与观测向量。、H：该生物系统的系统矩阵与观测矩阵。V、W：系统噪声向量和观测噪声向量。：系统真实状态与最优估计状态之间的偏差。:初始状态的最优平滑值。4 问题分析与求解4.1 问题一的求解该问题的实质是运用给定的观测数据，对系统参数进行辨识。由于该方程组为非线性微分形式，应用现有的线性系统参数辨识方法难以有效进行。而且由于传统差分法精度太差，难以满足高精度参数辨识的要求，需要寻求另外的有效方案。分析所给方程组，我们发现由所给微分方程组模型消去可得
9、 (1)分离变量后易得 (2)对上式进行积分得到： (3) (4)式（4）减式（3）可消去积分常数C，进而得到 (5)将DATA1所有数据代入（4）式，可得到关于的齐次线性方程组。 （j=0,.,4） (6)对于实际问题来说，该方程组的系数矩阵非奇异的概率几乎为1，从而使得该齐次线性方程组仅有零解。而问题1已经给出的具体数值，则可将（6）式化为如下形式的非齐次线性方程组：() (7) 式（7）为关于待定参数的超定方程组，因此不存在精确解，只存在某种意义下的最优解。本文采用最小二乘法求使得所有方程的误差平方和最小的参数解，即最小方差解。为此令 (8) (9) (10)由此可将式（7）写成如下形式
10、的矩阵方程 (11)设式（11）中的估计值为，则最小二乘估计的指标是使得各次观测值与由参数估计值确定的量测估计值之差的平方和最小，即 (12)而要使上式达到最小，需满足 (13)若系数矩阵H列满秩（实际问题中该条件始终是满足的），即正定，则的最小二乘估计为 (14)运行程序problem1.m即可得到待定参数（）的最小方差估计值。另一方面由所给模型可以看出、，并且题给DATA1.txt数据准确无误，因此待定参数的具体数值如表1所示。表1 根据DATA1及的数值确定的参数（）的值-2.0000060793290412.00021335922440-1.0000213305844810604.2
11、问题二的求解该问题与问题（1）的唯一差别是为未知，要求运用最少的数据对参数进行辨识。同上题一样，只需第一组数据即可确定参数及，分别为： (15)然而，由于未知，所以不能像上题一样把它转换为线性非齐次方程组进行求解，但是依然可以通过变换得到 (16)令 （） (17)其中为常数，则有： (18)只需将前4组数据带入上式即可得到一个线性非齐次的方程组（同样对实际问题其系数矩阵非奇异的概率为1），由此运行程序problem2.m容易确定出、的值分别为： (19)则原微分方程组可化为： (20)由于和（）只相差一个比例系数，所以只需要求出即可得出参数的值。下面通过微分方程的数值解法确定出生物种群数量与
12、待定参数之间的关系式。对于微分方程组的数值解法可以采用诸如梯形公式、辛普生公式、龙格库塔方法等经典方法来进行，分析表明，龙格库塔方法的截断误差可以通过积分步长进行较好的控制，当采用4级4阶龙格库塔公式时，数值积分方法的截断误差为，即的高阶无穷小。因而若取积分步长为时，其截断误差达到，可以较好的满足高精度参数估计的需要。另一方面，在先行64位微机上采用MATLAB软件进行计算时，可以不考虑其舍入误差。4级4阶龙格库塔法的计算公式如下： (21)其中： (22)然而当方程组中含有未知参数，采用龙格库塔方法求解微分方程组的数值解将难以进行具体的迭代运算，因而本文采用搜索寻优的方法来求解待定参数的具体
13、数值。具体方法是，首先在较大的范围内取较大的步长确定出的大概范围，然后逐步缩小搜索范围并减小步长来确定出满足精度要求的值。将搜索范围设定在之间，搜索步长取为1，运行程序，可以首先得到参数的粗略值为：在精确搜索时，根据粗略搜索的结果可将搜索范围设定为，搜索步长缩小为0.01，运行程序，可以进一步得到参数的值为：如此进行下去，不断缩小搜索范围和搜索步长，最后可确定出精度更高的的值，本文将的误差设定在之内，再次运行程序后，计算得到的值为：将满足精度要求的常数的值代入式（17）并利用式（19）结果，即可得到待定参数（）的值。在本问题中我们用到了四组数据完成了对参数的求解，结果如表2所示：表2 根据DA
14、TA1数值确定的参数（）值-14.164999189655101.4164873894973584.99175847016191-7.0826783142661410604.3 问题三的求解为了能够从含有误差的数据中估计出参数的最优值，采用控制系统中应用最为广泛的Kalman滤波方法进行。将非线性系统参数的辨识问题转换为各状态参数变量的最优状态估计，从而求得Kalman滤波意义下(其本质为最小二乘意义下)的参数最优辨识解具体处理方法如下。4.3.1 总体思想 此问题要求我们在我们在观测数据有误差的情况下，确定参数在某种意义下的最优解。 在有m组观测资料的情况下，我们的目标是求解最佳参数，使得观
15、测数据和方程组的解尽可能的接近，因此我们采用最小二乘思想建立非线性回归分析模型：其中 通过求解该最优化问题，确定出参数向量,在此参数中我们微分方程求解，然后得到实测数据和数值解法的残差平方和。4.3.1.1 算法流程图图2 模式搜索过程由于轴向搜索和模式搜索不断的调整的值，计算相应的，向着目标函数值减小的方向逐步前进，该算法能保证收敛，但计算量比较大，所以有时候是局部收敛。4.3.1.2 求解非线性常微分方程组一般可以看成是以下初值问题的解：在微分方程组解得存在唯一性条件保证后，设时间步长为又因为有其中均为二维向量。4.3.1.3 模式搜索法求解无约束优化1. 算法的基本思想算法从初始节点出发
16、，交替实施两种搜索：轴向搜索和模式搜索。轴向搜索依次沿n个坐标轴方向进行，用来确定新的基点和有利于函数值下降的方法。模式搜索沿着相邻的两个基点的连线方向进行，试图使函数值下降更快。2. 算法的步骤描述Step 1:初始化 给定初始点，初始步长，加速因子，缩减率，精度。Step 2:轴向搜索 如果，则令，转Step 3； 如果，则令，转Step 3； 否则，令。Step 3: 若，则令，转Step 2； 如果，转Step 4； 否则转Step 4.Step 4：模式搜索 令； 令转Step 2。Step 5： 如果，停止，得到点； 否则，令； 令转Step 2。该算法将轴向搜索和模式搜索的固定步
17、长改为用一维搜索的确定步长，算法仍然收敛。3. 算法图解图3 算法求解过程示意4.3.2 模型求解与数据分析根据以上理论，我们利用MatLab编制了该模型求解的程序。对于所给DATA2数据，得到Kalman滤波辨识过程中的，变化结果和各参数的协方差变化过程，分别如图 、图 所示。图中结果表明，滤波得到的，是渐进收敛的。改变系统给定的初值，滤波的稳定值没有变化。而且随着滤波次数的不断增加，各参数的协方差组成的协方差阵趋近于一个常值正定阵，且协方差均很小，进一步说明该滤波算法的稳定性和有效性。通过对程序中协方差矩阵COV的进一步分析，发现其中COV中有关、的协方差系数均为量级，说明对于、的估计比较
18、精确。图 4 利用DATA2得到的参数估计（辨识）过程图 5利用DATA3各辨识参数的协方差变化对于所给DATA3数据，同样按照上述方法，得到辨识过程中的，变化结果和各参数的协方差变化过程，结果分别如图 、图 所示。图 6利用DATA3得到的参数估计（辨识）过程图 7利用DATA3各辨识参数的协方差变化分析所得结果，我们发现该组中，参数辨识，趋于收敛，而，则出现周期性振荡现象，且通过各参数协方差变化过程的分析，认为、不能满足高精度参数辨识的要求。这是由于DATA3中采样周期约为DATA2的，采样点比较密集，观测资料有误差，采样数据易失效。因此辨识出的结果不尽理想。但通过对图4，图5，图6，图7
19、的综合分析，仍然可以看出，趋向收敛于统一值，表明该算法具有正确性和适用性。4.3.3 模型的改进受上述两组辨识结果的启发，以下通过模型改进，可以获得较好的辨识结果。具体方法为：在DATA3中，采用等间隔（间距为10）抽样的方法，将数据分为组，分别对每组利用上述方法进行参数辨识，然后将各组得到的参数进行加权平均，得到最终的辨识结果。该方法能够剔除辨识结果的随即误差，相当于对各组Kalman滤波意义下辨识出的最优的系统参数再次进行平滑，消除其随机误差。编程得到各组的参数辨识过程结果比较图和各辨识参数的协方差变化图，分别如图8，图9所示。图 8各组的参数辨识过程结果比较图图 9 各辨识参数的协方差变
20、化过程模型参数的校验：根据辨识出的，利用龙格库塔法及观测初值，进行数值积分，得到种群数目的标称值，与用扩展Kalman滤波算法得到的种群数目的最优估计值进行比较，如图10所示。图 10 龙格库塔方法与最优估计值的比较根据得到的精确的，采用固定点最优平滑算法，可以得出、的最优估计值，最终辨识出系统参数(利用DATA3)如表3所示。表3 利用数据DATA3得出的系统参数值1.87789685903883-0.0982064879995812-8.768090016198310.89518396879414712.8190546726550073.41087060535932利用DATA2得出的结果
21、表4如示。表4利用数据DATA2得出的系统参数值1.82454930199629-0.09306342741035-9.187905391459520.9189525683286512.8090546826550073.310890655359324.4 问题（4）的求解由于DATA4所给数据中对参数估计结果影响非常重要的时间变量也含有误差，所以无法直接沿用问题（3）中的求解方法来处理DATA4所提供的观测数据，这里如果将时间变量与种群数目的观测值、待定参数等同看待，均作为该食饵捕食者生物系统的状态变量，利用输出方程对时间变量进行观测，将能够借助扩展的Kalman滤波方法估计出各个状态变量（包
22、括时间变量的值），得到（）及时间变量在最小方差意义下的最优估计值。这样就可以利用固定点最优平滑算法来确定、在最小方差意义下的最优估计值。4.4.1 系统模型的建立与问题（3）的求解过程类似，所不同的是增加了时间状态变量，具体处理方法如下：设、，列出该生态系统非线性形式的状态方程： (23)令，表达式的含义与问题（3）相同，则式（46）可写为如下形式： (24)为系统误差（这里可假定为高斯白噪声）。由所给数据选择物种的数量及时间为观测量，即、，则该生物系统的观测方程即为： (25)式中为测量误差（这里仍然假定为高斯白噪声）。令，将式（24）和式（25）写成如下的矩阵形式： (26) (27)由于
23、方程（25）是非线性的，因此首先对（25）式线性化得到偏微分矩阵： （28）然后按照问题（3）改进算法进行Kalman滤波即可得出参数（）的估计值及物种的数量、和时间的滤波值，再采用固定点最优平滑得出、的最优估计值。4.4.2 模型的求解与数据分析根据4.4.1 节建立的数学模型，我们利用MatLab编制了该模型求解的程序。对于所给DATA4数据，得到Kalman滤波辨识过程中的，变化结果和各参数的协方差变化过程，分别如图11、图12所示。图11 分组滤波得到的参数辨识过程图 12分组滤波得到的参数协方差变化情况图中结果表明，滤波得到的，是渐进收敛的，而且随着滤波次数的不断增加，各参数的协方差
24、组成的协方差阵趋近于一个常值正定阵，且协方差均很小，说明问题三中提出的改进算法是很有效性的。与问题三中的处理方法相同，可以根据得到的待定参数，的精确值，采用固定点最优平滑算法，求出、的最优估计值，最终辨识出系统参数如表5所示。表5利用数据DATA4得出的系统参数值1.88391412547936-0.0978804717345565-8.815649634382370.89728529808289213.205360610418670.82496656898685 模型精度分析及评价5.1模型精度分析：对于问题一，我们利用全部六组无误差观测数据，根据极小范数最小二乘解求解，在已知时简单的比例运
25、算得到模型参数。此过程误差只能来自矩阵求逆过程和矩阵乘法运算的截断误差，提高计算过程中数据精度就能够提高模型参数的精度，直接来自真实的观测数据精度得到自然保证；对于问题二，利用四组数据得到，其精度依赖于数据运算中的数据精度。在求解的过程中，我们采用了在微分方程数值解中搜索观测数据从而确定时间间隔的方法。其精度受到微分方程数值解中步长和搜索算法的影响，减小步长，提高搜索深度均能够提高的精度，从而提高的精度，直接来自真实的观测数据精度得到自然保证；对于问题三，利用多组误差观测数据根据极小范数最小二乘得到，其精度依赖于数据运算中的数据精度和观测数据中的误差，的精度受到微分方程数值解中步长和搜索算法的
26、影响。在模型的改进中，我们增加了对于初值和常数的优化过程，由优化过程的最小均方误差目标，我们知道优化过程可以提高参数的准确度和精度；对于问题四，当观测时间存在误差时，的精度依赖于数据运算中的数据精度和观测数据中的误差。在时间间隔搜索过程，我们采用了预估周期校正算法，其数据精度可以通过搜索深度和步长进一步提高。以上我们对各个问题中模型参数的精度进行了定性分析，进一步我们可以通过函数的泰勒级数展开和误差传递进行定量分析。5.2模型评价如上分析该模型最大的优点是保证了参数估计的高精度，该模型的缺点是，建立的生态系统参数高精度估计模型，对于小幅度噪声具有一定的适应性，但是当噪声强度较大时，该模型得到的
27、结果精度显著下降。6模型的推广6.1生态过程周期的多解性在第二，三，四问的求解过程中，我们均采用了固定初始点搜索相邻点的搜索算法，在算法中我们均假设相邻两个数据来自真实生态过程的一个周期，此种假设下我们根据观测数据得到的生态过程周期是所有可能周期的最大值。这里我们给出生态过程周期的一般表达式。设该生态系统的真实演变过程周期为，而获取的数据其实是对真实演变过程的采样。这里就涉及了数据处理中一般的采样问题，根据奈奎斯特采样定理，当采样频率大于信号中最高频率时，可以根据采样数据恢复原始数据，进而从采样数据的周期得到信号的真实周期。但一般的数据观测过程不一定满足采样定理，因此同样的观测数据对应着多组真
28、实的周期。对于原始演变过程不同采样频率的采样可能得到相同的采样数据，这是一个一对多的关系，因此不能从采样数据得到准确的演变周期，只能得到最大演变周期，前提是每两个相邻数据均属于一个演变周期，另外可以得到演变周期的一般表达式。如下两组数据均是对于同一过程的观测，初始时刻相同，观测间隔不同。6.2 第三、四问模型优化过程中初始值和系统常数的循环优化过程在第三问模型优化改进中，我们提出了首先利用观测数据选取初值优化系统常数使得总的均方差达到极小，然后利用优化的系统常数进一步优化初始参数使得总的均方差达到新的极小。这一过程可以继续进行，即反复优化初始参数和系统常数，最后当两次优化结果的差异小于设定的高
29、精度无限时，停止循环过程，此时得到了优化和的最优结果。由于每一次优化过程均要按照龙格-库塔方法求解微分方程组的数值解，故时间复杂度较大。参考文献1吕学富，飞行器飞行力学，西北工业大学出版社，1994年10月2姜启源，数学模型（第二版），高等教育出版社，1992年4月3徐仲，张凯院，陆全，冷国伟，矩阵论简明教程，科学出版社，2005年6月4张肇炽，徐仲，代数与几何基础，高等教育出版社，2001年6月5王沫然，MATLAB与科学计算（第二版），电子工业出版社，2003年9月6崔祜涛 ，史雪岩 ，崔平远 ，栾恩杰，绕飞慢自旋小天体的航天器运动分析，航空学报，第25卷第1期，2004年1月7封建湖，车
30、刚明，聂玉峰，数值分析原理，科学出版社，2001年9月8王慧忠，张新全，Lotka-Volterra数学模型在草地管理中的应用，草业学报，第15卷，第1期，2006年2月附表1.最小二乘法求解A=zeros(5,3);A(1,1)=log(0.1374480266382216/60);A(1,2)=-log(11.750840650304518 /10);A(1,3)=10-11.750840650304518 ;A(2,1)=log(7.108705996120129/60);A(2,2)=-log( 3.4133367257849176/10);A(2,3)=10-3.4133367257
31、849176;A(3,1)=log( 04251595082899424/60);A(3,2)=-log( 20.80921881438798/10);A(3,3)=10-20.80921881438798 ;A(4,1)=log( 0.7182384931413827/60);A(4,2)=-log( 5.231982481945481 /10);A(4,3)=10- 5.231982481945481 ;A(5,1)=log(26.706603701525214/60);A(5,2)=-log(26.925221604818102/10);A(5,3)=10-26.925221604818
32、102;b=-0.1374480266382216-60. ,7.108705996120129-60.,0.4251595082899424-60., 0.7182384931413827-60, 26.706603701525214-60;vpa(inv(B*B)*B*b)；结果a = -2.0000135599512575028000000000000 12.000417823214720828400000000000 -1.00004236950021052016000000000002.图形的画法ezplot(-2*log(y)+0.2*y-12*log(x)+x+13.82,1,3
33、0,0.01,70); title(图4.1.2);text(26.925221604818102,26.706603701525214,rp)clearezplot(-2*log(y)+0.2*y-12*log(x)+x+13.82,1,30,0.01,70);title(图4.1.2);text(26.925221604818102,26.706603701525214,*)text(10,60,*)text(11.750840650304518,0.1374480266382216,*)text(3.4133367257849176,7.108705996120129,*) text(20.80921881438798,0.4251595082899424,*) text(5.231982481945481,0.7182384931413827,*)