基于改進序貫均勻設計的常微分方程組參數估計

2021-10-30 05:08:22王福昌賀財寶

濱州學院學報 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)表示按照得到的參數估計結果繪制的擬合曲線與觀測值的對比，可以看到估計出的參數符合觀測結果。