基于改進(jìn)序貫均勻設(shè)計(jì)的常微分方程組參數(shù)估計(jì)

王福昌,賀財(cái)寶

(防災(zāi)科技學(xué)院 基礎(chǔ)部，河北 廊坊 065201)

0 引言

均勻設(shè)計(jì)是我國(guó)數(shù)學(xué)家方開泰和王元提出的一種試驗(yàn)設(shè)計(jì)方法，特別適用于多水平多因素的試驗(yàn)設(shè)計(jì)復(fù)雜情形，主要通過均勻設(shè)計(jì)表進(jìn)行試驗(yàn)設(shè)計(jì)[1]。該方法自1978年被提出以來(lái)，引起了數(shù)學(xué)和統(tǒng)計(jì)學(xué)專家廣泛的關(guān)注[1-6]。

常微分方程組是一種應(yīng)用廣泛的數(shù)學(xué)模型，它的參數(shù)一般通過專家的專業(yè)知識(shí)確定。研究者一般假定模型參數(shù)已知時(shí)，討論方程的性質(zhì)，對(duì)于其反問題的研究相對(duì)較少。在實(shí)際應(yīng)用中，要想根據(jù)試驗(yàn)觀測(cè)數(shù)據(jù)篩選合適的模型，對(duì)它的參數(shù)進(jìn)行精確估計(jì)變得十分重要，有一些文獻(xiàn)對(duì)這一問題開展了研究[7-10]。針對(duì)最常見的獵物-捕食者模型的參數(shù)估計(jì)，人們提出了最小二乘逼近(LSA)法[7]、函數(shù)型數(shù)據(jù)分析法[8]和混合加速粒子群算法[9]。對(duì)于傳染病SIR模型，一般在最小二乘法準(zhǔn)則下用單點(diǎn)迭代法[10]。

以上傳統(tǒng)的迭代方法，有的對(duì)初值依賴,有的過于復(fù)雜，本文先把參數(shù)估計(jì)問題化為多變量?jī)?yōu)化問題，再用并行和精英保留策略改進(jìn)序貫均勻設(shè)計(jì)來(lái)估計(jì)微分方程組的參數(shù)，最后通過三個(gè)典型的應(yīng)用案例檢驗(yàn)方法的正確性和有效性。

1 基本假設(shè)和模型

假定常微分方程初值問題模型為

(1)

其中，y(t)=[y1(t),y2(t),…,yp(t)]T為關(guān)于時(shí)間的待求函數(shù)向量，θ∈Rs為模型參數(shù)，F(xiàn)(t,y(t),θ)為確定的非線性向量函數(shù)。再假設(shè)當(dāng)給定初值條件y(t0)=y0后，方程(1)存在唯一精確解。

設(shè)在時(shí)刻t1,t2,…,tn處y(t)的觀測(cè)向量值為y(t1),y(t2),…,y(tn)，預(yù)測(cè)向量值為h(t1,θ,y0),h(t2,θ,y0),…,h(tn,θ,y0)對(duì)應(yīng)的最小二乘估計(jì)

(2)

通過優(yōu)化模型(2)，即可得到參數(shù)θ,y0的最小二乘估計(jì)。

最小一乘估計(jì)在遇到離群值時(shí)，穩(wěn)健性更好，

也可給出M-估計(jì)形式[11]

ρ(·)可以取為Huber函數(shù)

其中c為常數(shù)。

有時(shí)優(yōu)化的目標(biāo)函數(shù)是預(yù)先設(shè)定的，如在連續(xù)時(shí)間確定性最優(yōu)控制問題中，優(yōu)化目標(biāo)是直接給出的，約束是一個(gè)常微分方程初值問題[12]

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

2 算法和應(yīng)用案例

均勻設(shè)計(jì)可以在參數(shù)候選解集內(nèi)均勻選點(diǎn)，利用較少的計(jì)算點(diǎn)來(lái)避免陷入局部最優(yōu)[1]。使用均勻設(shè)計(jì)方法需要先確定均勻設(shè)計(jì)表，而均勻設(shè)計(jì)表可以從網(wǎng)上下載(http://www.math.hkbu.edu.hk/UniformDesign/)或根據(jù)數(shù)學(xué)原理編程自動(dòng)生成[13]。

2.1 算法設(shè)計(jì)

設(shè)χ=[a,b]為Rs中一個(gè)超矩形，其中

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

即ai≤xi≤bi，i=1,2,…,s并設(shè)f(x)為χ上的連續(xù)函數(shù)，希望在試驗(yàn)區(qū)χ上找到x*，使得

則稱M為函數(shù)f(x)在區(qū)域χ上的全局最優(yōu)值，x*為χ上的全局最優(yōu)點(diǎn)。

這是一個(gè)典型的優(yōu)化問題，通過均勻設(shè)計(jì)來(lái)尋找全局最大點(diǎn)是一種可行的方法，其主要思路如下：給定試驗(yàn)區(qū)域χ上的一個(gè)設(shè)計(jì)

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

的一個(gè)設(shè)計(jì)點(diǎn)?？梢宰C明，在前面假設(shè)下，當(dāng)n→∞時(shí)，Mn→M，但這種方法的收斂速度還是很慢。Fang等提出的序貫均勻設(shè)計(jì)(Sequential Uniform Design Optimization，SNTO)[14]方法可以提高收斂速度。

設(shè)

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

為Cs=[0,1]上的設(shè)計(jì)，并設(shè)

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

2.2 改進(jìn)的序貫均勻設(shè)計(jì)方法

下面給出最小化目標(biāo)函數(shù)的序貫均勻設(shè)計(jì)算法。

步驟1 初始化。設(shè)t=0，χ(0)=χ,a(0)=a,b(0)=b。

步驟2 產(chǎn)生均勻設(shè)計(jì)。在試驗(yàn)區(qū)域χ(t)=[a(t),b(t)]尋找一個(gè)試驗(yàn)次數(shù)為nt的均勻設(shè)計(jì)P(t)，可以從網(wǎng)站下載，或者編寫程序自動(dòng)生成[13]。

步驟3 并行計(jì)算目標(biāo)函數(shù)新的近似值。由于每個(gè)點(diǎn)并沒有先后關(guān)系，因而可以用計(jì)算機(jī)的多個(gè)核同時(shí)計(jì)算目標(biāo)函數(shù)值。選取每個(gè)均勻設(shè)計(jì)點(diǎn)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 終止準(zhǔn)則。 令c(t)=(b(t)-a(t))/2，若maxc(t)<δ，其中δ為事先設(shè)定的很小的數(shù)，則χ(t)足夠小，且x(t)和Mt是可以接受的，此時(shí)終止算法，否則轉(zhuǎn)步驟5。

步驟5 更新試驗(yàn)區(qū)域并保留上次迭代中的最優(yōu)點(diǎn)。新的試驗(yàn)區(qū)域

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

其中

式中,γ稱為壓縮比。本次迭代中得到的最優(yōu)點(diǎn)進(jìn)入下次迭代，稱為精英保留策略。記t=t+1，轉(zhuǎn)步驟2。

2.3 應(yīng)用案例

例1 SIR模型參數(shù)估計(jì)。用英國(guó)傳染病監(jiān)測(cè)中心1978年發(fā)布的有關(guān)流感病人的統(tǒng)計(jì)數(shù)據(jù)來(lái)估計(jì)SIR模型的參數(shù)。該數(shù)據(jù)統(tǒng)計(jì)了兩周內(nèi)每天某一所男孩寄宿學(xué)校流感爆發(fā)和流行情況，總共有763個(gè)學(xué)生，從發(fā)生一個(gè)染病者開始統(tǒng)計(jì)染病學(xué)生的人數(shù)，總共統(tǒng)計(jì)了15天，每一天統(tǒng)計(jì)得到的染病人數(shù)分別為1、3、7、25、72、222、282、256、233、189、123、70、25、11、4[10]。由于該數(shù)據(jù)收集時(shí)間很短，學(xué)生沒有因病死亡，即人口總數(shù)是常數(shù),b=0,d=0,模型為

由假設(shè)，通過等價(jià)變形可以化為一個(gè)關(guān)于感染者的二階微分方程模型

這里只有兩個(gè)參數(shù)θ=(β,γ)，可以通過可視化序貫均勻設(shè)計(jì)方法來(lái)直觀了解參數(shù)尋優(yōu)過程，見圖1。根據(jù)文獻(xiàn)信息，設(shè)置參數(shù)的初始搜索區(qū)間為β∈[0.001,0.003]，γ=[0.4,0.5]，即設(shè)a=[0.001,0.4]，b=[0.003,0.5]。這里選用了100水平2因素的均勻設(shè)計(jì)表，均勻設(shè)計(jì)表通過Matlab自動(dòng)生成[12]。壓縮系數(shù)取為γ=0.9，誤差取為10-6，后面采用同樣設(shè)置通過多次嘗試篩選，迭代出來(lái)的參數(shù)估計(jì)為θ*=(0.0021,0.4820)，圖1(a)表示預(yù)測(cè)值與觀測(cè)值的殘差平凡和對(duì)著迭代次數(shù)增加而變化的情況，圖1(b)表示按照得到的參數(shù)估計(jì)結(jié)果繪制的擬合曲線與觀測(cè)值的對(duì)比，可以看到估計(jì)出的參數(shù)符合觀測(cè)結(jié)果。

(a)殘差平方和隨著迭代次數(shù)變化

(b)感染者隨時(shí)間變化分布與擬合曲線

注意觀察可以發(fā)現(xiàn)，圖1(a)的目標(biāo)函數(shù)值一開始并沒有隨著迭代次數(shù)嚴(yán)格遞減，這是因?yàn)闉榱吮苊庠缡飕F(xiàn)象，即過早陷入局部極小值，在設(shè)計(jì)程序時(shí)，迭代的前5步?jīng)]有采用精英保留策略。后面兩個(gè)例采用了相同的前5步不用精英保留策略設(shè)置。

例2 酶積累問題(enzyme effusion problem)的數(shù)學(xué)模型可描述為[16]

其中

p=[p1,p2,p3,p4],y0=[y1(0),y2(0)]

為待估參數(shù)，觀測(cè)數(shù)據(jù)見表1。

表1 酶積累問題的部分觀測(cè)數(shù)據(jù)

根據(jù)表1數(shù)據(jù)和公式(2),可構(gòu)造一個(gè)含6個(gè)自變量的目標(biāo)函數(shù)，優(yōu)化參數(shù)范圍取為

0.1≤p1≤0.5,2.5≤p2≤3,0.3≤p3≤0.5,

0≤u4≤0.15,21≤y1(0)≤23,30≤y2(0)≤38。

這里選用了100水平6因素的均勻設(shè)計(jì)表，均勻設(shè)計(jì)表通過Matlab自動(dòng)生成[13]。目標(biāo)函數(shù)隨著迭代過程而變化的情況見圖2(a)，可見很快接近收斂值。由序貫均勻設(shè)計(jì)得到的參數(shù)估計(jì)為

p1=0.3081，p2=2.6985，p3=0.4099，p4=0.1114，y1(0)=21.9886，y2(0)=34.5574，

目標(biāo)函數(shù)值4631，圖2(b)給出了根據(jù)估計(jì)參數(shù)用ode45()計(jì)算出的y1(t)的擬合曲線和數(shù)據(jù)分析。該問題在文獻(xiàn)[16]中獲得的最好結(jié)果為

p1=0.3190，p2=2.7010，p3=0.4190，p4=0.1031，y1(0)=22，y2(0)=39，

對(duì)應(yīng)的最小二乘殘差平方和為5076.6，大于這里的4631。

(a)目標(biāo)函數(shù)值隨迭代次數(shù)的變化

(b)由估計(jì)參數(shù)計(jì)算出的曲線和數(shù)據(jù)分布

例3一個(gè)7變量恒溫連續(xù)攪拌槽反應(yīng)器問題(continuously stirred tank reactor(CSTR) problem)的數(shù)學(xué)模型可描述為[17]

其中q=u4+u1+u2,u1,u2,u3,u4為方程參數(shù)。本例優(yōu)化函數(shù)比較特殊，是求

的最大值，該問題實(shí)質(zhì)上是一個(gè)最優(yōu)控制問題[12]。這里目標(biāo)函數(shù)PI需要求數(shù)值積分，可以用Simpson復(fù)化積分公式[15]。為便于和文獻(xiàn)比較，初值取為

x0=[0.1883,0.2507,0.0467,0.0889,0.1804,0.1394,0.1046][17]。

參數(shù)范圍取為

0≤u1≤20,0≤u2≤4,0≤u3≤6,0≤u4≤20。

這里選用100水平4因素的均勻設(shè)計(jì)表，均勻設(shè)計(jì)表通過Matlab自動(dòng)生成[13]。編寫程序運(yùn)行后可得

u1=12.0485，u2=3.8651，u3=0.8317，u4=11.8667，

得到PI最大值19.8408。而在文獻(xiàn)[16]中，

u1=11.4550，u2=4.5222，u3=0.6865，

設(shè)u4=6.0，對(duì)應(yīng)PI最大值19.0437，本文計(jì)算結(jié)果有一定的改進(jìn)。本題是求最大值，編程時(shí)把目標(biāo)函數(shù)加負(fù)號(hào)，變?yōu)槿樽钚≈怠D3(a)顯示的是最大值目標(biāo)函數(shù)的相反數(shù)隨著迭代次數(shù)增加而變化的情況，圖3(b)顯示了根據(jù)參數(shù)估計(jì)繪制的x1(t)～x7(t)數(shù)值解曲線。

(a)目標(biāo)函數(shù)隨迭代次數(shù)變化

(b)根據(jù)估計(jì)參數(shù)繪制變量的曲線

與單點(diǎn)迭代和搜索算法相比，序貫均勻設(shè)計(jì)可以從多個(gè)均勻設(shè)計(jì)點(diǎn)同時(shí)開始迭代和搜索，具有本質(zhì)上的并行性，因而效率會(huì)大大提高。

3 結(jié)論

針對(duì)含有非線性微分方程組初值問題的復(fù)雜參數(shù)估計(jì)問題，由于它們的目標(biāo)函數(shù)的梯度很難計(jì)算，因而常使用無(wú)須梯度的直接搜索算法。本文使用搜索性能較好且計(jì)算機(jī)實(shí)現(xiàn)比較方便的序貫均勻設(shè)計(jì)方法對(duì)這一難題進(jìn)行了探討，通過三個(gè)典型應(yīng)用，說明了方法的有效性。