研究生数学建模竞赛优秀论文B.doc

1、一、问题的重述考虑航天器在仅受到地球万有引力、航天器自身发动机作用力的作用下作平面运动，将地球和航天器视为质点，建立航天器运动的数学模型。显然这样的数学模型在精度上是远远不能满足实际需要的，在其他要求精确制导等有关高科技的实际问题中，我们都面临着类似的问题：我们必须建立高精度的数学模型，必须高精度地估计模型中的大批参数，因为只有这样的数学模型才能解决实际问题，而不会出现差之毫厘，结果却失之千里的情况。由于航天器的问题太复杂，本题仅考虑较简单的确定高精度参数问题。假设有一个生态系统，其中含有两种生物，即： A生物和B生物，其中A生物是捕食者，B生物是被捕食者。假设时刻捕食者A的数目为，被捕食者B

2、数目为，它们之间满足以下变化规律：初始条件为：其中为模型的待定参数。通过对此生态系统的观测，可以得到相关的观测数据。要利用有关数据，解决以下问题：1) 在观测数据无误差的情况下，若已知，求其它5个参数？2)若也未知，至少需要多少组观测数据，才能确定参数？3) 在观测资料有误差（时间变量不含有误差）的情况下，确定参数 在某种意义下的最优解，并与仿真结果比较，进而改进数学模型。4) 假设连观测资料的时间变量也含有误差，确定参数在某种意义下的最优解。二、航天器运动模型的建立考虑航天器在仅受到地球万有引力、航天器自身发动机作用力的作用下作平面运动，将地球和航天器视为质点，由理论力学可知，一个刚体在空间

3、的运动可以看作质心的移动,因此可以应用质心运动定理来研究刚体质心的移动规律。以地球中心为原点，建立直角坐标系，航天器绕地球飞行，可以出现在该直角坐标系中四个象限的任意一个之内。平面直角坐标系如图1。符号说明如下:航天器在x方向的速度航天器在y方向的速度万有引力，航天器发动机作用力，为控制变量万有引力与x轴正方向的夹角航天器发动机作用力与x轴正方向的夹角初始时刻，航天器初始位置航天器x方向初速度航天器y方向初速度航天器受的万有引力方向指向地球中心（原点），航天器受推力的方向与x轴正方向成角。将和投影到该直角坐标系上，见图1图1 航天器受力分解图其中， 初始条件为，都是关于时间的位置函数，航天器在

4、x方向的分速度即对时间t求导，航天器在y方向的分速度即对时间t求导，航天器在x方向的加速度即对时间t求二阶导，航天器在y方向的加速度即对时间t求二阶导，根据牛顿第二定律有方程（3）和（4）。由此建立的航天器模型如下：显然这样的数学模型在精度上是远远不能满足实际需要的，在其他要求精确制导等有关高科技的实际问题中，我们都面临着类似的问题：我们必须建立高精度的数学模型，必须高精度地估计模型中的大批参数，因为只有这样的数学模型才能解决实际问题，而不会出现差之毫厘，结果却失之千里的情况。这时所建立数学模型的精度就成了数学模型的生命线。例如上述问题中的航天器还要受到地球质量分布不均匀所引起的摄动力，大气阻

5、力，日、月及其它星球的摄动引力的影响，以及航天器发动机为调整航天器自身姿态运作时作用力的影响。这样不但数学模型十分复杂，而且在这些数学模型中还要涉及到许多重要的参数，如地球的引力场模型就有许多待定参数。不仅如此，在对航天器进行测量时，还涉及到观测站的地理位置以及设备的系统误差等参数。为此人们要设法利用长期积累的丰富的观测资料，高精度确定这些重要的参数。由于航天器的问题太复杂，下面本题仅考虑较简单的确定高精度参数问题。三、捕食者与被捕食者生态系统问题的分析题中假设有一个生态系统，含有两种生物，A生物和B生物， A生物是捕食者，B生物是被捕食者。假设时刻捕食者A的数目为，被捕食者B数目为，它们之间

6、满足以下变化规律：初始条件为：该模型中，捕食者独自存在时死亡率，；被捕食者对捕食者的供养能力；是被捕食者的独立生存增长率，；是捕食者掠取被捕食者的能力，。2这个方程就是生态系统中被捕食者与捕食者的volterra模型，为模型的待定参数。对于该模型理论上不存在解析解，因此我们不能通过参数拟合确定模型的参数。Volterra模型在给定参数和初始值的情形下可以采用数值积分获得任意时间点的数值解。根据volterra模型进行一些公式推导如下：两个方程相除得：移项得：两边积分：得到相轨方程： 移项得：该式右边为只与系统初始状态有关，令易知，将（5）式代入（4）式得到 式中方程两边同除以C，得： 在（7）

7、式中，令，得到 这一方程体现了Volterra模型中两个变量之间的变化关系，我们称此方程为相轨方程。进一步研究相轨方程，可以发现Volterra模型中两个变量呈现周期性变化。第一问，对来说，相轨方程是一个4未知数的方程， DATA1中有6组数据，用6组数据确定4个可以采用极小范数最小二乘解。又因为已知，C可求，从而可求，由于各观测值真实准确，可取DATA1中任意一组数据，不失一般性，我们取第一组数据为初始值。第二问，我们可以证明参数C与系统周期成反比，由参考资料可以知道volterra模型中的含义，从而确定C的正负性，在未知的情况下，求C可以从C与时间的关系入手，我们先在DATA1的6组数据中

8、取4组算出，然后设计一个的新生态系统，以无误差的观测数据DATA1为准，设计搜索算法找到与DATA1中x，y值极为接近的数值点，找到对应的观测时间，得到观测间隔，这个观测间隔与DATA1中已知的观测间隔一起可以求出C，从而得到，同第一问。第三问，用所有数据求得极小范数最小二乘解，可以确定。经过与第二问类似的方法获得C。进一步求出 ，可取DATA2中观测初始时间的值，这一套就是我们所求的最小二乘意义下的最优解。将这一组参数带入volterra模型，获得各观测点上的仿真结果。通过与观测结果比较，我们发现误差普遍较大。于是我们改进了参数估计模型，改为求取均方误差意义下的最优解。获得了较好的效果。第四

9、问，我们采取了第三问中的改良算法，以求取使x,y,t三者均方误差最小的参数组为目标，进行计算，然而误差较大。经过判断，我们认为这是由于时间和数据均存在误差导致搜索结果不够精确，我们改进搜索算法，结果大为改善。四、模型的建立及求解问题一，已知，求其它5个参数结合DATA1.TXT中6组无误差的观测数据（包括了观测时刻、A生物数目、B生物数目 ），（7）式含有4个未知数，而题中提供了6组数据，写为矩阵形式即：记为矩阵形式，其中，矛盾方程组的惟一极小范数最小二乘解为，采用极小范数最小二乘解，得到的值，如表1所示表1 最小二乘的的值0.14471560254457-0.01447146313960-0

10、.868317924024250.07236037259773因为，由题意，而从volterra模型本身出发，是捕食者掠取被捕食者的能力，所以利用DATA1中数据算出的，所以C0，这一问中 ，把C代入，得到表2 已知， 的值对于，来说，由观测数据，已知，而未知，所以，可以是观测数据中任意一组，的值。不失一般性我们取，由DATA1知道，所以，。问题二，未知，至少需要几组数据才能确定的值1. 的求解由问题一的分析可知，为了求取，至少需要四个方程，即四组观测数据。列为线性方程可以解得，2. 观测系统时间变换分析对于同样系统的观测，当选取的观测时间起点和观测时间单位不同时，得到同一物理系统的两种不同时

11、间观测结果。不失一般性，我们将两种观测之间的时间变换关系表示为，其中表示两种观测系统的时间单位数量，反映了两个观测系统时间起点上的差异，而反映了两个观测系统观测时间单位的比例，线性关系是由时间的均匀性确定的。定理一：对于参数完全相同的生态系统，其系统常数与观测时间单位长度存在反比例关系，即。证明：在volterra模型中：把，代入得： （8）构建对于上述生态过程的两个观测系统，其时间轴变换关系为，为常数。（8）即为 （9）对（9）进行推导 （10）所以，为常数，结合，得到，这就证明了在参数完全相同的生态系统，其系统常数与观测时间间隔参数存在反比例关系，即。在volterra模型中，对的推导与上

12、面相同，结论也相同。3求解在确定的情况下，我们只要得到系统常数就可以确定生态系统参数。对于一个生态系统，当确定时，由方程可知相轨线是完全确定的。对于观测数据中的相邻两组数据和，其演变过程遵循系统方程即选择作为起始点根据系统方程演化到。当值不同时，从演化到的过程不同。下面我们在定理一中证明值和演化时间存在反比关系。我们首先用DATA1中的4组数据确定表3 的值0.144715602544570.01447146313960 0.868317924024250.072360372597731，模型的建立与求解龙格-库塔方法求解微分方程对于Volterra模型，没有显式的符号解，因此我们采用四阶龙格

13、库塔方法求解常微分方程组的数值解。求解方法介绍如下：volterra模型可写为 （11）令下标表示步数，则解此方程组的欧拉方法为 （12）引进向量记号，则式（11）与式（12）可分别写成此时，常用的四阶龙格库塔方法取形式采用龙格库塔方法，我们可以求出微分方程的数值解，数值解十分密集，在图2上表现为连续曲线。图2 龙格-库塔数值解与DATA1数据对照2，模型的建立我们的目标是使用最少数目的观测数据组，即目标： 约束：为DATA1中无误差观测间隔，为已知，为新系统观测间隔，这个新观测间隔由后面介绍的搜索算法可以获得，所以C可求，由DATA1中任意4组数据解出，由C已求出，则可求，验证的精度是用求出

14、的再建立新模型v_new，采用DATA1观测间隔得到观测值，q为该观测值与DATA1观测值的相对误差，这在结果验证中详细介绍。3，搜索算法我们采用一种搜索算法寻找与DATA1中无误差数据极为接近的高精度，以此确定新的观测间隔。由于观测数据无误差，我们可以选择DATA1中任意一组数据作为微分方程的初值，根据四阶龙格-库塔算法求解得到微分方程的数值解。不失一般性，我们选取第一组数据作为初值，即。以DATA1中的为起点进行龙格-库塔数值解的搜索。 搜索算法：第一步：预估生态系统的周期；第二步：在预估周期附近搜索真实周期；第二步：以为步长，生成1000组数据，搜索其中与目标点距离最近的点；在此点前后两

15、点构成的区间内重复上述搜索。随搜索深度上升，获得点的精度随之上升，我们在计算中考虑到精度和效率，选择搜索深度为三，即精确到 4，问题的求解为了计算，首先在龙格-库塔数值解基础上搜索出第一组x，y值，得到观测时间，然后继续搜索第二组x，y值，得到对应的观测时间，这时两观测时间间隔可求，如表4，表4 新系统观测间隔9.788112399999599.787980399999489.788139999999559.788111599999599.78790879999960记表4中各个观测间隔为，它为后一观测时间减去前一观测时间的值。在已知数据DATA1中，各个观测值与观测间隔如表5所示，表5 DA

16、TA1中无误差的观测数据及观测间隔 none10.60.0.111.7508406503045180.13744802663822160.13.41333672578491767.1087059961201290.1000000000000000420.809218814387980.42515950828994240.099999999999999965.2319824819454810.71823849314138270.126.92522160481810226.706603701525214记表5中各个观测间隔为，它为后一观测时间减去前一观测时间的值。由定理一，我们从理论上知道，结合表

17、4，表5，又因为，结合表2，的值即可确定，因为我们的值与无误差观测值有极小差距，我们的用表示。表6 (即)的值14.164836461043451.4164741404531484.991398119686477.08266991315792对于，来说，由观测数据，已知，而未知，所以，可以是观测数据中任意一组，的值。我们取，由DATA1知道，所以，。由此，我们通过构建一个新的观测系统找到了C与观测间隔t的对应关系，并且除了确定需要的DATA1中的4组数据之外，没有再需要DATA1中的其它数据，因此，至少需要4组数据就可以确定的值。5，结果验证为了验证我们的结果，对于得到的表6中的的值，加上，可

18、以再重新构建一个volterra模型（模型v_new），将DATA1中的观测时间代入模型v_new，得到，表7 模型v_new中取与DATA1相同观测间隔对应的x，y值 NAN10600.111.750933046976590.137451352070610.13.413376988757657.108234559668320.1000000000000000420.809288951375380.425167798549810.099999999999999965.232025537817290.718176598277300.126.9252566485123026.70655846347

19、881定义相对误差，其中x，y为DATA1中无误差的观测数据，得到相对误差表如表8所示。 表8 ，值与DATA1的相对误差表（单位：）0.049972095034430.011604495058510.078629840047040.241941078854410.11795781068136-0.663181811240250.033704767117530.194991754885340.08229360851517-0.861759216174480.01301519248942-0.01693889904760由此可见我们得到的，值具有5个数量级以上的精度，它与无误差观测值是极为接近的

20、。所以我们得到的的值是高精度的，如果需要我们可以加大搜索深度，提高数据精度。问题三，在观测资料有误差（时间变量不含有误差）的情况下，确定参数在某种意义下的最优解，并与仿真结果比较，进而改进模型。1，误差数据预处理：当观测数据有误差时，我们可以进行如下数据预处理以部分去除数据中误差的干扰。S1：根据全体数据对，按照极小范数最小二乘求解生态系统参数S2：因为只要四组数据就可以求解方程组，从而得到一组系统参数。但是由于误差的影响，使得这些参数与总体参数的偏差程度不同。我们可以设定偏差门限，从而去除含有较大误差的数据组。S3：对于剩余的数据对，选择一组作为初始值采用前面提到的搜索算法搜索相邻的下一组值

21、，得到观测时间间隔；S4：计算上一步计算得到的观测时间间隔的中位数，根据预先设定比例门限，去除时间间隔不属于中位数某一定义的邻域的数据对。经过以上步骤，在一定程度上去除了含有较大误差的数据组，减少了数据量，对于后续的优化过程具有重要作用。2，模型的建立及求解（1）DATA2系统由第二问求解过程，我们知道的求解只与数据点有关。我们可以利用第一步数据预处理后的全体数据方程做极小范数最小二乘解，得到。由定理一参数C与系统观测间隔成反比，而DATA2系统有给定的观测间隔，所以与第二问类似，建立的虚拟系统。假设DATA2中观测时间相邻的两点在同一周期内（第一点为初始值），则让新系统采用与DATA2系统相

22、同的初始值，经过与第二问类似的龙格-库塔计算及搜索，搜索的目标是x，y在新旧系统中对应值的均方差最小，DATA2系统的相邻两点的第二点反映到新系统内就得到定理一中所示的新系统对应于DATA2系统的观测时间，从而得到新系统观测间隔，即：表9 新系统与DATA2系统的对应关系DATA2系统新系统第一个观测点 观测间隔观测值x，y观测间隔搜索值第二个观测点根据定理一，经过DATA2中150组数据的运算，C可求（），对这150组数据采用最小二乘解，可求。DATA2中150组数据点和龙格-库塔数值解的相轨图如图3，图3 DATA2中150组数据点和DATA2系统的龙格-库塔数值解我们的目标是让x，y的均

23、方误差最小，由此建立目标规划模型目标： 约束：通过这个模型，我们可以直接求出满足它的C值，这个值是没有经过优化的，DATA2系统可直接解出它的极小范数最小二乘解：，根据，最优解如表7所示，（我们将记为），我们定DATA2中的初始点表10 DATA2系统的最优解， （2）DATA3系统DATA3中对建立与之相同的目标规划模型，也直接求出满足它的C值，这个值是没有经过优化的，DATA3系统也可直接解出它的极小范数最小二乘解：， ， 根据，最优解如表8所示，（我们将记为），我们同样定DATA3中的初始点表11 DATA3系统的最优解3，所求模型与仿真结果比较DATA2中，对求出的最优建立生态系统，把

24、无误差的观测时间代入该系统，可以发现在无误差观测时间下x，y的值与DATA2中x，y的值有较大误差，目标7.37696260317394，同样，对DATA3，对求出的最优建立生态系统，把无误差的观测时间代入该系统，可以发现在无误差观测时间下x，y的值与DATA3中x，y的值的均方差33.47612623110470，也非常大，所以表7，表8中的值是非常不精确的，下面设计改良模型的方法。4，模型的改良（1）DATA3系统首先将DATA3中的初始点固定，修正新系统的C，取区间，采用二分法进行搜索，当二分法的搜索精度保证C精确到时就停止搜索，找到满足该搜索目标的作为改良的C值，这时的均方差1.192

25、06591140184， 确定值，反过来再修正，这里就包含了迭代的思想，由于C一定，且不变，则一定，我们认为的修正值在DATA3的初始值附近，给定的（同C）为区间，对于各个不同的 与构成不同系统，寻找使最小的系统，这通过对调用MATLAB语句做最小二乘曲线拟合来实现，这个系统对应的初值就是的修正值。 ，1.18427007187202综上可见：未优化的均方差33.47612623110470首先优化C值，优化后1.19206591140184， 然后优化初值，优化后 ，1.18427007187202P值在第一次优化之后的显著下降，在第二次优化之后进一步下降，从而可以知道这种优化方法是非常好的

26、。至此，我们可以建立改良模型如下：目标： 约束：（2）DATA2系统DATA2系统采用与DATA3系统相同的优化方法，计算得 ，0.25392540697850问题四，假设连观测资料的时间变量也含有误差，确定参数在某种意义下的最优解。当观测资料的时间变量也含有误差时，使得观测时间的间隔含有更大的误差。如上所述，求解时，相当于求解相轨线方程，这与时间变量没有关系，因此这四个参数的求解精度仅受观测数据误差的影响。同问题三的求解过程，我们在得到时，选取的新的观测系统，则这是一个除初值外均确定的系统。设为搜索得到的时间间隔，该时间间隔与这两组数据有关，搜索结果为时间间隔与之间的对应关系，如果时间没有误

27、差，我们可以直接计算，但是由于时间误差的存在，直接计算将会把时间上的误差以倒数形式传递到上，当较小时即使很小的时间误差都会对结果造成巨大的影响。因此我们对以上常数计算方法进行改进，提出了一下新的计算方法以降低误差的影响。新的常数计算方法为。即选取连续若干组数据，每次均是从前一组搜索相邻下一组数据，从上式可以看出经过这样的处理后，分母上是一段比较长的时间间隔，两个时间的误差对间隔的影响将会大大降低，从而使得计算得到的常数更加稳定。实际处理中还可以多次分组，对已经比较稳定的常数进行最小二乘拟合以进一步降低误差水平，提高数据精度。通过以上方法确定了常数，下面考虑优化初值，方法同问题三求解过程中的循环

28、优化方法。因为时间，观测值均含有误差，因此我们进行参数优化高精度求解的目标不应该仅仅包含观测数据的均方误差，还应该包含时间的均方误差，即均方误差和为，为DATA4系统时间观测值，为观测时间的估计值。根据以上分析，我们建立时间含有误差下的单目标非线性规划模型目标： 约束：其中为DATA4系统时间观测值，为新系统时间观测值时间含有误差情形下时间估计算法：S1：对于给定的数据组DATA4，进行数据预处理；S2：利用预处理结果数据计算；S3：预估数据周期；S4：搜索时间间隔，并估计；S5：对和初值进行优化；S5：代入参数和初始值计算均方误差；DATA4系统可直接解出它的极小范数最小二乘解：（，）通过这

29、个模型，我们可以直接求出满足它的C值。经过优化获得因为计算时间问题，我们没有进行初始值优化。最终以第一组参数为初始值获得的仿真数据中，两变量和时间的均方差分别为Pxy=0.54933666827330，P=0.00711664554174综合为P=0.77691196210573该值为取第一组参数为初始值条件下最优。五、模型精度分析及评价模型精度分析：对于问题一，我们利用全部六组无误差观测数据，根据极小范数最小二乘解求解，在已知时简单的比例运算得到模型参数。此过程误差只能来自矩阵求逆过程和矩阵乘法运算的截断误差，提高计算过程中数据精度就能够提高模型参数的精度，直接来自真实的观测数据精度得到自然

30、保证；对于问题二，利用四组数据得到，其精度依赖于数据运算中的数据精度。在求解的过程中，我们采用了在微分方程数值解中搜索观测数据从而确定时间间隔的方法。其精度受到微分方程数值解中步长和搜索算法的影响，减小步长，提高搜索深度均能够提高的精度，从而提高的精度，直接来自真实的观测数据精度得到自然保证；对于问题三，利用多组误差观测数据根据极小范数最小二乘得到，其精度依赖于数据运算中的数据精度和观测数据中的误差，的精度受到微分方程数值解中步长和搜索算法的影响。在模型的改进中，我们增加了对于初值和常数的优化过程，由优化过程的最小均方误差目标，我们知道优化过程可以提高参数的准确度和精度；对于问题四，当观测时间

31、存在误差时，的精度依赖于数据运算中的数据精度和观测数据中的误差。在时间间隔搜索过程，我们采用了预估周期校正算法，其数据精度可以通过搜索深度和步长进一步提高。以上我们对各个问题中模型参数的精度进行了定性分析，进一步我们可以通过函数的泰勒级数展开和误差传递进行定量分析。模型评价如上分析该模型最大的优点是保证了参数估计的高精度，该模型的缺点是，建立的生态系统参数高精度估计模型，对于小幅度噪声具有一定的适应性，但是当噪声强度较大时，该模型得到的结果精度显著下降。六、模型的推广一 生态过程周期的多解性。在第二，三，四问的求解过程中，我们均采用了固定初始点搜索相邻点的搜索算法，在算法中我们均假设相邻两个数

32、据来自真实生态过程的一个周期，此种假设下我们根据观测数据得到的生态过程周期是所有可能周期的最大值。这里我们给出生态过程周期的一般表达式。设该生态系统的真实演变过程周期为，而获取的数据其实是对真实演变过程的采样。这里就涉及了数据处理中一般的采样问题，根据奈奎斯特采样定理，当采样频率大于信号中最高频率时，可以根据采样数据恢复原始数据，进而从采样数据的周期得到信号的真实周期。但一般的数据观测过程不一定满足采样定理，因此同样的观测数据对应着多组真实的周期。对于原始演变过程不同采样频率的采样可能得到相同的采样数据，这是一个一对多的关系，因此不能从采样数据得到准确的演变周期，只能得到最大演变周期，前提是每

33、两个相邻数据均属于一个演变周期，另外可以得到演变周期的一般表达式。如下两组数据均是对于同一过程的观测，初始时刻相同，观测间隔不同。对于以上两个观测，当时，两组观测结果完全相同，其中为整数。根据观测过程我们有如下关系式：，其中为两种观测情形下得到的观测数据周期，为生态系统最大演变周期，为生态系统的一般周期表达式。二 第三、四问模型优化过程中初始值和系统常数的循环优化过程在第三问模型优化改进中，我们提出了首先利用观测数据选取初值优化系统常数使得总的均方差达到极小，然后利用优化的系统常数进一步优化初始参数使得总的均方差达到新的极小。这一过程可以继续进行，即反复优化初始参数和系统常数，最后当两次优化结果的差异小于设定的高精度门限时，停止循环过程，此时得到了优化和的最优结果。由于每一次优化过程均要按照龙格-库塔方法求解微分方程组的数值解，故时间复杂度较大。22