确定Lotka-Volterra生态系统模型高精度参数的研究(精)

第 37 卷第 14 期 2007 年 7 月 数学的实践与认识

MA THEMA T ICS IN PRACT ICE AND THEORY Vo l137No 114 July, 2007

确定 Lotka -Volterra 生态系统模型高精度参数的研究 唐静波 , 胡智渊 , 张彦琼 指导教师 : 任力锋

(中南大学生物医学工程研究院 , 湖南长沙 410078

摘要 : 研究确定 Lo tka2Vo lterra 生态系统模型的高精度参数估计问题 . 利用周期性 , 先对测量数据进行预

处理 ; 然后用三种不同的方法构造了误差函数 , 进行非线性最小二乘法参数估计 ; 再用计算机仿真对其进行

验证 . 结果表明该方法能够有效地解决高精度参数估计中消除测量数据误差的问题 .

关键词 : 非线性最小二乘法 ; 参数估计 ; 周期性 ; 平均值法 1问题提出

收稿日期 : 2007203201

在 “神舟六号” 载人航天宇宙飞船、 人造地球卫星等航天器围绕地球在轨运行的过程中 , 必须建立高精度的数学模型和高精度地估计模型中的大批参数 . 因为只有这样的数学模型 才能解决实际问题 , 而不会出现差之毫厘 , 结果却失之千里的情况 . 为此人们要设法利用 长

期积累的丰富的观测资料 , 高精度确定这些重要的参数 .

由于航天器的问题太复杂 , 下面仅考虑较简单的确定高精度参数问题 . 假设有一个生态系统 , 其中含有两种生物 , 即 : A 生物和 B 生物 , 其中 A 生物是捕食者 , B 生物是被捕食者 . 假设 t 时刻捕食者 A 的数目为 x ( t , 被捕食者 B 数目为 y ( t , 它们 之间

满足以下变化规律 : x ′ ( t = x ( t [a 1 + a 2y ( t ] y ′ ( t = y ( t [a 3 + a 4x ( t ] 初始条件为 : x ( t 0 = a 5 y ( t 0 = a 6 (1

其中 a k (1 F k F 6 为模型的待定参数 .

该模型为著名的 Lo tka2Vo lterra 方程 [ 1— 2 ]. 通过对此生态系统的观测 , 可以得到相关的 观测数据 . 为了说明问题 , 本文仅以 DA TA 1. TXT , DA TA 4. TXT 为例 [ 3 ].

问题 1. 1在观测数据无误差的情况下 , 至少需要多少组观测数据 , 才能确定参数 a k (1 F k F 6 ? 有关数据见数据文件 :DA TA 1. TXT

问题 1. 2在观测数据 x 、 y 、 t 均有误差的情况下 , 试建立数学模型 , 确定参数 a k (1 F k F 6 在某种意义下的最优解 , 并与仿真结果比较 . 有关数据见数据文件 : DA TA 4. TXT 2模型求解

Lo tka2Vo lterra 模型的假设 :

1 假设是一个包括捕食与被捕食两个种群的生态系统 , 捕食者靠被捕食者而生存 , 系 统与外界没有种群交换关系 .

2 假设种群自身的增殖变化速率与 t 时刻种群的数量成比率 (马尔萨斯定律 . 2. 1方程通解法

对于问题 1. 1, 可由 Lo tka2Vo lterra 方程求得方程的通解 , 然后利用其通解 , 构建误差函 数 E

1, 由于参数 a k 越精确 , E

1 的值越接近 0. 由此可以转化为对函数 M in E 1 的非线性最小值

无约束优化模型 . 对于此非线性最小值无约束优化模型的求解 , 我们采用非线性最小二乘 法

对参数 a 1~a 4 和参数 C 进行估计 [ 4— 5 ]. 另外根据测量资料 DA TA 1. TXT 得知 x ( t 0 = a 5 = 10, y ( t 0 = a 6 = 60.

Lo tka2Vo lterra 方程通解 : a 1 ln y - a 3 ln x + a 2y - a 4x = C (2

误差函数 : E

1 = (a 1 ln y - a 3 ln x + a 2y - a 4x - C 2 (3 M inE 1 = Σ N i = 0

[a 1 ln y i ( t - a 3 ln x i ( t + a 2y i ( t - a 4x i ( t - C ]2 (4 MA TLAB 实现过程 :

1 用 Lo tka2Vo lterra 方程通解 , 构建误差函数 E 1 (3 ,M inE 1 (4.

2 用非线性最小二乘法对参数 a 1~a 4 和参数 C 进行估计 (M at lab fm in search 函数 .

3 用估计得到的参数对 Lo tka2Vo lterra 方程进行求解 (M at lab ode45 函数 . 4 将仿真结果与原始数据进行比较 .

仿真结果与原始数据非常吻合 , 见图 1. 图中横坐标表示捕食者 x ( t , 纵坐标表示被捕 食者 y ( t ; 原始数据用红色‘ 3 ’表示 , 仿真数据用蓝色实线表示 .

我们做了一个统计 , 分别使用 5 组、 4 组、 3 组数据进行参数估计的结果见表 1. 结果表 示 : 至少需要 4 组观测数据可以确定参数 , 但计算结果不够精确 , 建议至少使用 5 组数据 . 图 1DA TA 1. TXT 数据与仿真数据图 2DA TA 4. TXT 数据与仿真数据

对比图 (通解法 对比图 (通解法

对于问题 1. 2, 由于测量资料 DA TA 4. TXT 的 x ( t , y ( t 和时间变量 t 均有误差 , 故 需要对

这些数据进行预处理 . 我们首先对原始数据进行插值 : t i = 0: 0. 01: 15, 求出 t i 时刻对应 的 x ( t ,

y ( t , 以消除时间变量 t 的误差 . 然后利用周期解定理 (见 2. 2 的定理 1 , 对插值后的 数据按周期

T 进行分段 , 共 9 段 (每段长 162 个数据点 , 共 1458 个数据点 . 将 1 至 9 段中对应时间点

的 x ( t ,

y ( t 加起来 , 求其平均值 x q i ( t 、 y q

i ( t , i = 1, 2, ⋯ , 162, 以减少数据 x ( t , y ( t 的误差 . 78 数 学 的 实 践 与 认 识 37 卷

表 1用 5 组、 4 组、 3 组数据对 Lo tka2Vo lterra 方程进行参数估计统计表

数据组 5 组 4 组 3 组 模 型 参 数

a 1 - 0. 23239389370315 - 0. 25718506931796 - 0. 60526159001990 a 2 0. 02323002030168 0. 02581146824179 0. 27384481371427 a 3 1. 39566452999218 1. 54218980966463 1. 45825956312549 a 4 - 0. 11628571769750 - 0. 12853240689509 - 0. 15877718474521 SS 残差平方和

1. 85807e- 007 2. 05723e- 007 M S= SSön 平均残差平方和

3. 7161400e- 008 5. 143075e- 008 仿真曲线和原始数据 偏差很大 , 故不推荐 只使用 3 组数据 . 利用误差函数 E

1 和函数 M inE

1, 采用非线性最小二乘法对参数 a 1~a 4 和参数 C 进行估

计 , 参数 a 5, a 6 可由测量资料 DA TA 4. TXT 中直接得到 . 仿真结果见图 2. x i ( t =

Σ 9 j = 1

x ( j - 1 ×162+ i ( t 9 , y i ( t = Σ 9 j = 1

y ( j - 1 ×162+ i ( t 9 ( i = 1, 2, ⋯ , 162 (5 2. 2平均值法

定理 1设 x ( t , y ( t 是 Lo tka2Vo lterra 方程 (1 的一个周期解 , 具有周期 T > 0. 定 义 x

和 y 的平均值为

x q = 1 T ∫ T x ( t dt , y q = 1 T ∫ T y ( t dt , (6 那么 , x q = - a 3 öa 4 及 y q = - a 1 öa

2, 换句话说 , x ( t 和 y ( t 的平均值就是平衡值 (图 3. 图 3Lo tka2Vo lterra 模型的轨迹

证明 用 x ( t 去除 (1 的第一个等式的两边 , 得到 x ′ x

= a 1 + a 2y , 于是 1 T ∫ T x ′ x d t = 1 T ∫ T

[a 1 + a 2y ]dt (7 令 ∫ T x ′ x

d t = ln x (T - ln x (0 , 因为 x (T = x (0 , 于是这个等式等于零 . 因此

14 期唐静波 , 等 : 确定 Lo tka2Vo lterra 生态系统模型高精度参数的研究 79 1

T ∫ T a 2y ( t dt = 1 T ∫ T

- a 1d t = - a 1, (8 于是 , y q = - a 1 öa

2. 类似地 , 用 y ( t 去除 (1 的第二个等式的两边并从 0 到 T 进行积分 , 可得 x q = - a 3 öa 4.

对于问题 1. 2, 由于测量资料 DA TA 4. TXT 的 x ( t , y ( t 和时间变量 t 均有误差 , 故 需

要对这些数据进行预处理 . 我们首先对原始数据进行插值 : t i = 0: 0. 01: 15, 求出 t i 时刻 对

应的 x ( t , y ( t , 以消除时间变量 t 的误差 . 然后利用周期解定理 (定理 1 , 对插值 后的数据

按周期 T 进行分段 (每段长 162 个数据点 , 共 9 段 , 共 1458 个数据点 , 并分别将 1 至 9 段内

的 x ( t , y ( t 的 162 个数据点加起来求平均值 x q i ( t , y q

i ( t , i = 1, 2, ⋯ 9, 以减少系统变量 x ( t , y ( t 的误差 .

利用周期解定理 (定理 1 中的种群的平均值模型 , 构造误差函数 E 2 (11 , 并采用非线性

最小二乘法对参数 a 1~a 4 进行估计 , 仿真结果见图 4.

x i ( t = Σ 162 j= 1 x ( i- 1 ×162+ j ( t 162 y i ( t = , Σ 162 j= 1 y ( t 162 ( i = 1, 2, ⋯, 9 (9 种群的平均值模型: ( i- 1 ×162+ j x q = - a3 öa 4 y q = - a1 öa 2 ] a3 + a4x q = 0 a1 + a2y q = 0 (10 误差函数: E 2 = [a1 + a2y q i ( t ]2 + [a3 + a4x q i ( t ]2 (11 M inE 2 = Σ N i= 1

{[a1 + a2y q i ( t ]2 + [a3 + a4x q i ( t ]2} (12 图4 图5 DA TA 4. TXT 数据与仿真 数据对比图(平均值法 DA TA 4. TXT 数据与仿真数据 对比图(四阶龙格2库塔法 四阶龙格2库塔法 2. 3 四阶龙格 库塔法 由于已知在DA TA 4. TXT 中测量数据x ( t , y ( t , t 均有误差, 故首先对原始数据 进行 插值, 求出ti 时刻( ti = 0: 0. 001: 15 对应的x ( t , y ( t , 在这里时间t 的增量, 即步长h 取 01001. 然后根据四阶龙格2库塔法公

式, 求出Lo tka2Vo lterra 方程中x ( t , y ( t 的 增量d x , 80 数 学 的 实 践 与 认 识37 卷 d y , 并利用d x , d y 构建误差函数E 3 (14 , 用非线性最小二乘法对参数a1～ a4 进行参数估 计, 以减少系统变量x ( t , y ( t 和时间t 的误差. 仿真结果见图5. d x n= x n+ 1- x n= h 6 (k1+ 2k2+ 2k 3+ k4 k 1= x n (a1+ a2y n k 2= x n+ h 2 a 1+ a 2 y n+ h 2 k1 k 3= x n+ h 2 a 1+ a 2 y n+ h 2 k2 k 4= (x n+ h [a1+ a2 (y n+ hk3 ] dy n= y n+ 1 - y n= h 6 (k1+ 2k2+ 2k 3+ k4 k1= y n (a3+ a4x n k 2= y n+ h 2 a 3+ a 4 x n+ h 2 k 1 k 3= y n+ h 2 a 3+ a 4 x n+ h 2 k 2 k4= (y n+ h [a3+ a4 (x n+ hk 3 ] (13

误差函数: E 3 = [d x n - d x 3 n ]2 + [d y n - d y 3 ]2 (14 M inE n 3 = Σ N i= 0 {[d x n - d x 3 n ]2 + [d y n - d y ]2} (15 d x n 3 n , d y n: 根据(13 计算得到增量 d x 3 n , d y : 对DA TA 4. TXT 数据进行插值, 然后计算得到的增量. 3 n 3 模型验证 我们对所有模型验证的方式都采用把原始数据与计算机仿真数据进行比较、计算残差 平方和及平均残差平方和的方法. 计算结果显示, 残差平方和都在可接受误差范围内. 4 模型评价和推广 本文根据Lo tka2Vo lterra 模型的周期性, 先对测量数据进行预处理以消除测量数据的 误差, 然后用三种不同的方法构造误差函数进行非线性最小二乘法参数估计. 计算机仿真 的 结果表明该方法能够有效地解决高精度参数估计中消除测量数据误差的问题. 参考文献: [ 1 ] [ 2 ] [ 3 ] [ 4 ] [ 5 ] W ILL IAM E LUCA S. 微分方程模型[M ]. 长沙: 国防科技大学出版社, 1998. 徐克学. 生物数学[M ]. 北京: 科学出版社, 2001. h ttp: ööwww. shumo. com 第三届全国研究生数学建模竞赛B 题. 刘庆上, 赵捍东, 王芳. 用优化方法求解最小二乘法问题[J ]. 弹箭与制导学报, 2004, 24 (3 : 83—86. 张志涌. 精通MA TLA TB6. 5 版[M ], 北京: 北京航空航天大学出版社, 2003. Study on Determ ination of High Accuracy Parameter in Lotka-Volterra Artif ic ial Ecosystem Model TAN G J ing2bo, HU Zh i2yuan, A dviso r: REN L i2feng ZHAN G Yan2qiong ( Inst itute of B iomedical Engineering of Cent ral South U niversity, Changsha Hunan 410078, Ch ina

14 期唐静波, 等: 确定Lo tka2Vo lterra 生态系统模型高精度参数的研究81 Abstract: Abstract To study the est imat ion p roblem on determ inat ion of h igh accuracy parameter in Lo tka2Vo lterra art ificial eco system model. A t first, the

measured datum w ere p ret reated by periodicity, and then const ruct th ree different erro r funct ion w h ile p rocessing nonlinear least square method parameter est imat ion; finally, to verify the results by computer simulat ion. In conclusion, it is effect ive to elim inate erro rs of the measured datum in h igh accuracy parameter est imat ion by th ismethod. Keywords: Keywords method nonlinear least squaremethod; parameter est imat ion; periodicity; average number 期刊简介 本刊主要刊登数学的最新的理论成果, 及其在工业、农业、环境 保护、军事、教育、科研、经济、金融、管理、决策等工程技术、自 然科学 和社会科学中的应用成果、方法和经验. 主要任务是沟通数学工作 者与其他科技工作者之间的联系, 推动应用数学在我国的发展, 为 四 化建设作贡献. 主要栏目: 数学建模、管理科学、工程、应用、问题研究、知识与进 展、学科介绍、方法介绍、高等数学园地、数学史、研究简报、书刊 评 介、简讯. 注: ① 投稿一式二份, 原稿自己保存, 编辑部不退还投稿的稿件. ② 编辑部不接收网上投稿. 82 数 学 的 实 践 与 认 识 37 卷__