基于改进序贯均匀设计的常微分方程组参数估计

2021-10-30王福昌贺财宝

滨州学院学报 2021年4期

王福昌,贺财宝

(防灾科技学院基础部，河北廊坊 065201)

0 引言

均匀设计是我国数学家方开泰和王元提出的一种试验设计方法，特别适用于多水平多因素的试验设计复杂情形，主要通过均匀设计表进行试验设计[1]。该方法自1978年被提出以来，引起了数学和统计学专家广泛的关注[1-6]。

常微分方程组是一种应用广泛的数学模型，它的参数一般通过专家的专业知识确定。研究者一般假定模型参数已知时，讨论方程的性质，对于其反问题的研究相对较少。在实际应用中，要想根据试验观测数据筛选合适的模型，对它的参数进行精确估计变得十分重要，有一些文献对这一问题开展了研究[7-10]。针对最常见的猎物-捕食者模型的参数估计，人们提出了最小二乘逼近(LSA)法[7]、函数型数据分析法[8]和混合加速粒子群算法[9]。对于传染病SIR模型，一般在最小二乘法准则下用单点迭代法[10]。

以上传统的迭代方法，有的对初值依赖,有的过于复杂，本文先把参数估计问题化为多变量优化问题，再用并行和精英保留策略改进序贯均匀设计来估计微分方程组的参数，最后通过三个典型的应用案例检验方法的正确性和有效性。

1 基本假设和模型

假定常微分方程初值问题模型为

(1)

其中，y(t)=[y1(t),y2(t),…,yp(t)]T为关于时间的待求函数向量，θ∈Rs为模型参数，F(t,y(t),θ)为确定的非线性向量函数。再假设当给定初值条件y(t0)=y0后，方程(1)存在唯一精确解。

设在时刻t1,t2,…,tn处y(t)的观测向量值为y(t1),y(t2),…,y(tn)，预测向量值为h(t1,θ,y0),h(t2,θ,y0),…,h(tn,θ,y0)对应的最小二乘估计

(2)

通过优化模型(2)，即可得到参数θ,y0的最小二乘估计。

最小一乘估计在遇到离群值时，稳健性更好，

也可给出M-估计形式[11]

ρ(·)可以取为Huber函数

其中c为常数。

有时优化的目标函数是预先设定的，如在连续时间确定性最优控制问题中，优化目标是直接给出的，约束是一个常微分方程初值问题[12]

x′=f(x,u),x=x(t)=xi(初值)。

2 算法和应用案例

均匀设计可以在参数候选解集内均匀选点，利用较少的计算点来避免陷入局部最优[1]。使用均匀设计方法需要先确定均匀设计表，而均匀设计表可以从网上下载(http://www.math.hkbu.edu.hk/UniformDesign/)或根据数学原理编程自动生成[13]。

2.1 算法设计

设χ=[a,b]为Rs中一个超矩形，其中

a=[a1,a2,…,as],b=[b1,b2,…,bs],

即ai≤xi≤bi，i=1,2,…,s并设f(x)为χ上的连续函数，希望在试验区χ上找到x*，使得

则称M为函数f(x)在区域χ上的全局最优值，x*为χ上的全局最优点。

这是一个典型的优化问题，通过均匀设计来寻找全局最大点是一种可行的方法，其主要思路如下：给定试验区域χ上的一个设计

P={xk,k=1,2,…,n}，

的一个设计点。可以证明，在前面假设下，当n→∞时，Mn→M，但这种方法的收敛速度还是很慢。Fang等提出的序贯均匀设计(Sequential Uniform Design Optimization，SNTO)[14]方法可以提高收敛速度。

设

P0={yk,k=1,2,…,n}

为Cs=[0,1]上的设计，并设

xki=ai+(bi-ai)yki,i=1,2,…，s，xk=(xk1,…,xks)，k=1,2,…,n,

2.2 改进的序贯均匀设计方法

下面给出最小化目标函数的序贯均匀设计算法。

步骤1 初始化。设t=0，χ(0)=χ,a(0)=a,b(0)=b。

步骤2 产生均匀设计。在试验区域χ(t)=[a(t),b(t)]寻找一个试验次数为nt的均匀设计P(t)，可以从网站下载，或者编写程序自动生成[13]。

步骤3 并行计算目标函数新的近似值。由于每个点并没有先后关系，因而可以用计算机的多个核同时计算目标函数值。选取每个均匀设计点x(t)∈P(t)∪{x(t-1)}，使得

Mt=f(x(t))≤f(y)，∀y∈P(t)∪{x(t-1)},

式中,x(-1)表示空集，x(t)和Mt分别为x*和f(x*)的近似值。

步骤4 终止准则。令c(t)=(b(t)-a(t))/2，若maxc(t)<δ，其中δ为事先设定的很小的数，则χ(t)足够小，且x(t)和Mt是可以接受的，此时终止算法，否则转步骤5。

步骤5 更新试验区域并保留上次迭代中的最优点。新的试验区域

χ(t+1)=[a(t+1),b(t+1)]，

其中

式中,γ称为压缩比。本次迭代中得到的最优点进入下次迭代，称为精英保留策略。记t=t+1，转步骤2。

2.3 应用案例

例1 SIR模型参数估计。用英国传染病监测中心1978年发布的有关流感病人的统计数据来估计SIR模型的参数。该数据统计了两周内每天某一所男孩寄宿学校流感爆发和流行情况，总共有763个学生，从发生一个染病者开始统计染病学生的人数，总共统计了15天，每一天统计得到的染病人数分别为1、3、7、25、72、222、282、256、233、189、123、70、25、11、4[10]。由于该数据收集时间很短，学生没有因病死亡，即人口总数是常数,b=0,d=0,模型为

由假设，通过等价变形可以化为一个关于感染者的二阶微分方程模型

这里只有两个参数θ=(β,γ)，可以通过可视化序贯均匀设计方法来直观了解参数寻优过程，见图1。根据文献信息，设置参数的初始搜索区间为β∈[0.001,0.003]，γ=[0.4,0.5]，即设a=[0.001,0.4]，b=[0.003,0.5]。这里选用了100水平2因素的均匀设计表，均匀设计表通过Matlab自动生成[12]。压缩系数取为γ=0.9，误差取为10-6，后面采用同样设置通过多次尝试筛选，迭代出来的参数估计为θ*=(0.0021,0.4820)，图1(a)表示预测值与观测值的残差平凡和对着迭代次数增加而变化的情况，图1(b)表示按照得到的参数估计结果绘制的拟合曲线与观测值的对比，可以看到估计出的参数符合观测结果。