琚新剛，廉飛宇，董 樂(lè)，張 元，葛宏義，蔣玉英

(1.河南財(cái)政金融學(xué)院 計(jì)算機(jī)與信息技術(shù)學(xué)院 計(jì)算機(jī)軟件與理論重點(diǎn)學(xué)科，河南 鄭州 450000；2.河南工業(yè)大學(xué) 糧食信息處理與控制教育部重點(diǎn)實(shí)驗(yàn)室，河南 鄭州 450001)

0 引 言

生物模擬系統(tǒng)的參數(shù)估計(jì)是一個(gè)NP hard問(wèn)題，具有多模態(tài)和病態(tài)的特征，被估計(jì)的參數(shù)通常有多個(gè)真實(shí)解[1]。梯度優(yōu)化算法，由于對(duì)初始值敏感度高，通常不能得到真實(shí)解[2]。啟發(fā)式算法，已被證明能夠解決NP hard問(wèn)題，能夠克服多峰性問(wèn)題，并得到全局最優(yōu)解[3]。

生物模擬系統(tǒng)的參數(shù)估計(jì)問(wèn)題，有無(wú)跡卡爾曼濾波(unscented Kalman filter，UKF)[4,5]、模擬退火(simulated annealing，SA)算法[6]、螢火蟲(chóng)算法[7]、遺傳算法[8]等。

而改進(jìn)的粒子群優(yōu)化算法(particle swarm optimization，PSO)，使用了一種分解技術(shù)，可以使算法更易于求得最優(yōu)解，可使粒子在全局解附近只有較小的擾動(dòng)，從而避免了粒子大的波動(dòng)和適應(yīng)度函數(shù)值的惡化。因此，與普通的PSO算法相比，本算法有著更低的均方根誤差(root mean squared error，RMSE)。

為了驗(yàn)證算法的有效性，將改進(jìn)的PSO算法應(yīng)用于一個(gè)仿真的生物模型的參數(shù)估計(jì)(該生物模型的數(shù)據(jù)來(lái)源于一個(gè)E-Coli系統(tǒng)的真實(shí)數(shù)據(jù))。比較改進(jìn)PSO算法與迭代無(wú)跡卡爾曼濾波器(iterated unscented Kalman filter，IUKF)、SA和PSO算法的RMSE值，結(jié)果表明，改進(jìn)的PSO算法具有更小的RMSE值，從而能夠更精確地估計(jì)生物系統(tǒng)仿真模型的參數(shù)。

1 數(shù)學(xué)建模

假定生物系統(tǒng)的非線性動(dòng)力學(xué)模型如下

(1)

x∈N是代謝物向量，u∈m是生物分子的濃度向量，即自變量，w∈q是參數(shù)向量。對(duì)該模型的求解,就是要找到最適合模型的未知參數(shù)的解。

用來(lái)描述生物現(xiàn)象的建模框架有多種可供使用。諸如Lotka-Volterra[9]、S-systems[10]等建?？蚣苤校x擇了S-system來(lái)表達(dá)式(1)所描述的非線性模型

(2)

式中：αi>0和βi>0為比例系數(shù)，g和h是活躍度。N是代表新陳代謝的因變量的數(shù)目，m是作為系統(tǒng)輸入的獨(dú)立變量的數(shù)目。其中，x=[x1,…,xN]T為狀態(tài)向量(即代謝物)，u=[xN+1,…,xN+m]T為輸入向量(即生物分子的濃度向量)，w=[α1,…,αN,β1,…,βN,g11,…,gN,N+m,h11,…,hN,N+m]T為參數(shù)向量，q=N(2+2(N+m)) 是未知參數(shù)的最大數(shù)目。

2 模型求解方法

有許多方法可用于估計(jì)以上非線性方程的未知參數(shù)。為了將本方法與其它啟發(fā)式方法進(jìn)行性能對(duì)比，選擇了模擬退火方法。此外，為了進(jìn)一步分析這種啟發(fā)式方法，還選擇了一種基于無(wú)跡卡爾曼濾波器的非啟發(fā)式方法。除了所提出的算法進(jìn)行描述，還將描述以上提到的這些對(duì)比方法。

2.1 迭代無(wú)跡卡爾曼濾波器(IUKF)

UKF是卡爾曼濾波器用于非線性建模的改進(jìn)框架[11]。在此濾波器中，過(guò)程和檢測(cè)方程可能是非線性的，因而需要非線性卡爾曼濾波器。UKF算法采用無(wú)跡變換作為非線性變換來(lái)計(jì)算隨機(jī)變量的統(tǒng)計(jì)量。

UKF算法如下：

式(1)的離散狀態(tài)空間模型

x[k+1]=F(x[k],u[k],w)

(3)

由定義y[k]=x[k+1],k=0,…,N-1， 上述方程可以寫(xiě)成如下的非線性過(guò)程方程

y[k]=F(x[k],u[k],w)

(4)

式中：x[k] 和u[k] 為輸入，y[k] 為輸出，w為待估計(jì)的具有維數(shù)q的未知參數(shù)。為了估計(jì)參數(shù)，定義如下?tīng)顟B(tài)空間

w[k+1]=w[k]+r[k]y[k]=F(x[k],u[k],w[k])+e[k]

(5)

其中，前者為過(guò)程方程，r為過(guò)程噪聲。后者則是由輸入和測(cè)量噪聲e[k] 驅(qū)動(dòng)的測(cè)量方程。

UKF算法的收斂速度較慢，濾波器可能不收斂于真實(shí)參數(shù)的主要原因，是不當(dāng)?shù)某跏贾捣峙?，?dǎo)致了卡爾曼增益的過(guò)早飽和。因此，較小的濾波器增益導(dǎo)致參數(shù)估計(jì)較小的變化，會(huì)使算法的收斂變得非常緩慢。

為解決這個(gè)問(wèn)題，開(kāi)發(fā)了IUKF算法。在本算法中，RMSE計(jì)算公式如下

(6)

2.2 模擬退火(SA)算法

SA算法是Scott Kirkpatrick等提出的一種元啟發(fā)式算法。該算法基于材料的加熱和冷卻原理，從一個(gè)先驗(yàn)解出發(fā)，通過(guò)過(guò)程的重復(fù)不斷改進(jìn)算法本身，直到得到問(wèn)題的最優(yōu)解。算法由一個(gè)內(nèi)循環(huán)和一個(gè)外循環(huán)組成。內(nèi)部循環(huán)使用本地搜索產(chǎn)生解空間中的最終解，并基于概率準(zhǔn)則更新該解。外循環(huán)持續(xù)地降低過(guò)程的溫度，這個(gè)溫度會(huì)影響內(nèi)循環(huán)的性能。該算法在求解高溫問(wèn)題的初始階段，具有較好的搜索解空間的能力，避免了非凸問(wèn)題的局部解。通過(guò)降低溫度，算法顯示了良好的探測(cè)能力。

在高溫時(shí)的第一次迭代中，與代價(jià)函數(shù)的非適當(dāng)值相反，內(nèi)循環(huán)中出現(xiàn)壞解的概率較高。這個(gè)特性可防止算法在局部響應(yīng)中收斂。隨著溫度的降低，在外環(huán)的最后迭代中，壞解的概率降低了，算法以較高的概率選擇了最合適的解。

2.3 粒子群算法優(yōu)化(PSO)

PSO算法是Kennedy和Eberhart[12]提出的一種基于種群的隨機(jī)優(yōu)化方法。PSO模擬了鳥(niǎo)類、昆蟲(chóng)群或魚(yú)群的活動(dòng)，在這些活動(dòng)中，個(gè)體被稱為粒子，群體被稱為粒子群。每個(gè)粒子都有一個(gè)位置P和一個(gè)速度V，并分散在搜索空間中。粒子在搜索空間中按照它們最知曉的位置分布，同時(shí)也按照整個(gè)群體最知曉的位置分布。迭代t中粒子i第j維的速度為

(7)

式中：W為慣性權(quán)重，通常在0.9～1.2范圍內(nèi)，c1和c2是設(shè)為2的常數(shù)，r1和r2為區(qū)間[0,1]內(nèi)的隨機(jī)參數(shù)。

第i個(gè)粒子的新位置更新公式如下

Pij(t)=Pij(t-1)+Vij(t)

(8)

PSO算法是一種快速簡(jiǎn)單的方法，適用于非凸NP hard問(wèn)題，如生物路徑建模。雖然它在搜索解空間上有很好的表現(xiàn)，但在最優(yōu)解的鄰域中搜索最優(yōu)解上存在不足。針對(duì)這一不足，通常采用一些啟發(fā)式算法或?qū)⑵渑c其它方法結(jié)合來(lái)解決。在粒子群算法中嵌入一種新的啟發(fā)式算法，提高了算法的尋優(yōu)能力。在該算法中，在一個(gè)序列平臺(tái)中使用規(guī)范化的PSO開(kāi)發(fā)了一種分解算法，該算法稱為基于分解的PSO算法(particle swarm optimization algorithm base on decomposition，PSOBD)，由兩部分組成。

第一部分，是一個(gè)標(biāo)準(zhǔn)的PSO，算法在最優(yōu)解的鄰域內(nèi)找到一個(gè)解，即全局最優(yōu)解。為了提高第一部分盡可能接近最優(yōu)解的能力，在PSO算法的每次迭代中，線性減小式(7)中的慣性權(quán)重，最終值為0.1。當(dāng)粒子在最佳解附近時(shí)，使用這種技術(shù)，可以在最后的迭代中，提供更精確的解。

第二部分，采用分解技術(shù)，偽代碼如下：

begin

-以全局最佳的方式重新定位所有粒子的位置

-在全局最佳位置的鄰域，在粒子位置的第j維上散開(kāi)粒子

-使用PSO算法求第j維的最優(yōu)值

-更新全局最佳位置

end

end

如上述偽代碼所述，在算法的第二部分，PSO初始化為隨機(jī)粒子，然后通過(guò)迭代，找到最優(yōu)解。原本在每一次迭代中，粒子通過(guò)兩個(gè)臨時(shí)機(jī)制來(lái)更新自身的值。現(xiàn)在不用整個(gè)種群，只用其中一部分最優(yōu)粒子的相鄰值，相鄰值的極值就是最優(yōu)值。內(nèi)部循環(huán)使用規(guī)范的PSO，試圖找到所有決策變量的最優(yōu)值，而將其它決策變量的值設(shè)為常數(shù)值。在第二部分，通過(guò)分解問(wèn)題到每個(gè)決策變量，把PSO算法從多維搜索空間集中到利用單維空間上。

當(dāng)其它決策變量作為固定的常數(shù)值時(shí)，PSOBD算法能夠在最佳解的附近，搜索到?jīng)Q策變量的最優(yōu)值。這種方法對(duì)所有決策變量是可重復(fù)的。在外部循環(huán)中，重復(fù)整個(gè)過(guò)程，以確保避開(kāi)陷入局部最優(yōu)。該方法可以產(chǎn)生更高精度的最優(yōu)解。

3 仿真與實(shí)驗(yàn)

3.1 仿真結(jié)果與分析

采用Matlab2016a進(jìn)行仿真。

通過(guò)兩個(gè)綜合仿真來(lái)說(shuō)明改進(jìn)的粒子群算法的性能。在第一個(gè)仿真中，假設(shè)沒(méi)有附加噪聲。在第二個(gè)仿真中，考慮了信噪比為20的加性高斯白噪聲。數(shù)據(jù)通過(guò)以下模型生成

(9)

參數(shù)如下

ω=[α1,α2,α3,α4,β1,β2,β3,β4,g13,g15,g21,g32,g41,h11,h22,h33,h34,h44]

(10)

表1給出了上述模型的參數(shù)，模型通過(guò)4個(gè)隨機(jī)初始值，生成4個(gè)合成的新陳代謝系統(tǒng)。這4種合成狀態(tài)的時(shí)間分布由式(9)得到，分布如圖1所示。4種合成代謝時(shí)間分布中圖1(a)和圖1(d)相似，但圖1(a)變化較緩，圖1(b)隨時(shí)間有升和降的交替變化，而圖1(c)有一個(gè)單峰的變化。根據(jù)式(9)，隨機(jī)產(chǎn)生4個(gè)初始值，產(chǎn)生4種合成代謝。這些數(shù)據(jù)是用于估計(jì)模型的參數(shù)，利用模型預(yù)測(cè)數(shù)據(jù)和真實(shí)數(shù)據(jù)計(jì)算RMSE。

表2～表5給出了IUKF、SA、PSO和提出的PSOBD算法的仿真結(jié)果。在這些表中，給出了每種算法被估參數(shù)的1000次仿真平均值和相應(yīng)的平均RMSE值。PSOBD算法在沒(méi)有附加噪聲時(shí)RMSE值為0.197，分別比IUKF算法、SA算法和PSO算法下降60.12%、71.16%和49.87%，而在有附加噪聲(SNR=20)時(shí)，RMSE值為0.304，分別比IUKF算法、SA算法和PSO算法下降46.95%、64.02%和38.83%?？梢钥闯觯琍SOBD算法對(duì)兩種仿真場(chǎng)景下的未知參數(shù)估計(jì)，都有較好的效果，該方法的RMSE遠(yuǎn)小于其它3種算法。在兩個(gè)仿真實(shí)驗(yàn)中，RMSE值降幅達(dá)39%～71%。仿真結(jié)果表明，與普通粒子群算法相比，PSOBD算法的參數(shù)估計(jì)能力更強(qiáng)，與IUKF算法和SA算法相比，PSOBD算法具有更精準(zhǔn)的參數(shù)估計(jì)能力。圖2給出了表2～表5中每種算法在有噪聲和無(wú)噪聲情況下1000次仿真的RMSE值變化情況。

表1 模擬S系統(tǒng)中選擇的參數(shù)

圖1 4種合成狀態(tài)的時(shí)間分布

相對(duì)于SA、IUKF和PSO方法，PSOBD算法在無(wú)噪聲和有噪聲場(chǎng)景下的平均百分比分別提升到60.06%和48.72%。

3.2 實(shí)驗(yàn)結(jié)果與分析

為了比較算法在實(shí)驗(yàn)場(chǎng)景中的性能，利用了Cad系統(tǒng)中的真實(shí)數(shù)據(jù)集。Cad系統(tǒng)是大腸桿菌條件應(yīng)力響應(yīng)模塊之一，僅在低pH和富含賴氨酸的環(huán)境下被激活。下面的S系統(tǒng)可以用來(lái)模擬這種現(xiàn)象

表2 IUKF算法兩次仿真的估計(jì)參數(shù)

表3 SA算法兩次仿真的估計(jì)參數(shù)

表4 PSO算法兩次仿真的估計(jì)參數(shù)

表5 PSOBD算法兩次仿真的估計(jì)參數(shù)

圖2 在無(wú)噪聲和有噪聲兩種情況下1000次仿真運(yùn)行的RMSE值

(11)

利用上述方法對(duì)系統(tǒng)參數(shù)進(jìn)行了估計(jì)，這些參數(shù)將用于由式(11)生成的Cad系統(tǒng)。為了比較算法的結(jié)果，采用了RMSE作為指標(biāo)衡量生成的時(shí)間分布圖與真實(shí)數(shù)據(jù)集的誤差。表6顯示了計(jì)算得到的RMSE。

由估計(jì)參數(shù)得到生成的時(shí)間分布和真實(shí)時(shí)間分布如圖3所示。

由圖3可見(jiàn)，將分解技術(shù)應(yīng)用于一般的粒子群算法中，

表6 4種算法在實(shí)際實(shí)驗(yàn)中的RMSE

使得PSOBD算法在全局解附近的搜索活動(dòng)減少，從而避免了較大的跳躍對(duì)結(jié)果的影響。因此，如前所述，與普通的PSO算法相比，PSOBD算法具有更好的逼近性能。在真實(shí)數(shù)據(jù)集的實(shí)驗(yàn)中，PSOBD算法的RMSE為0.807，相對(duì)于SA、IUKF和PSO方法，分別改進(jìn)了20.34%、27.10%和11.42%，與IUKF、SA和PSO算法的比較結(jié)果說(shuō)明了該方法的良好性能。

綜合仿真和真實(shí)數(shù)據(jù)集上的實(shí)驗(yàn)，測(cè)試結(jié)果表明，本算法的均方根誤差比其它對(duì)比方法分別平均減少了55.16%和19.62%。

圖3 由估計(jì)參數(shù)得到的時(shí)間分布和真實(shí)時(shí)間分布

4 討 論

討論了4種不同的啟發(fā)式算法在具有大量參數(shù)的非線性生物系統(tǒng)中的參數(shù)估計(jì)問(wèn)題，啟發(fā)式算法在參數(shù)估計(jì)中的應(yīng)用很廣泛。IUKF可提供卡爾曼濾波器數(shù)學(xué)意義上的解。除了IUKF外，還考慮了兩種元啟發(fā)式方法：SA算法和PSO算法。它們是求解NP hard問(wèn)題最優(yōu)解的兩種強(qiáng)有力的方法，而PSOBD算法具有更好的尋優(yōu)性能。實(shí)驗(yàn)結(jié)果表明，該算法在全局解附近的跳變較小，相對(duì)于其它3種算法的RMSE更小。從RMSE的角度看，PSOBD算法優(yōu)于其它3種方法。

在仿真中，建立了一個(gè)帶有17個(gè)參數(shù)的模型，并在兩種不同信噪比的情況下分析了4種算法的參數(shù)估計(jì)性能，并用估計(jì)的和真實(shí)的參數(shù)計(jì)算RMSE。PSOBD的RMSE比其它3種算法都要小。

在實(shí)驗(yàn)中，利用計(jì)算機(jī)輔助設(shè)計(jì)系統(tǒng)對(duì)一個(gè)真實(shí)的生物代謝系統(tǒng)進(jìn)行了建模。該模型與真實(shí)場(chǎng)景相似，我們對(duì)模型參數(shù)應(yīng)用以上的4種算法進(jìn)行了估計(jì)。實(shí)驗(yàn)結(jié)果表明，對(duì)于IUKF、SA、PSO和PSOBD算法在CPU為i5 4590，內(nèi)存8 G條件下，每個(gè)算法的運(yùn)行時(shí)間分別為343 s、305 s、287 s和416 s，PSOBD算法運(yùn)行時(shí)間較長(zhǎng)，表明PSOBD方法有進(jìn)一步改進(jìn)的空間。

5 結(jié)束語(yǔ)

提出了一種基于改進(jìn)粒子群算法的非線性生物模型的未知參數(shù)估計(jì)算法PSOBD。通過(guò)采用分解技術(shù)，把PSO算法從多維空間集中到單維空間上，當(dāng)把其它決策變量作為常數(shù)時(shí)，PSOBD算法能夠搜索到?jīng)Q策變量的最優(yōu)值。該算法對(duì)所有決策變量是可重復(fù)的，在外部循環(huán)中，可以避開(kāi)陷入局部最優(yōu)，產(chǎn)生更高精度的最優(yōu)解。將PSOBD算法得到的結(jié)果與IUKF算法、SA算法和一般的粒子群優(yōu)化算法進(jìn)行了比較，結(jié)果表明，與其它3種方法相比，PSOBD算法在速度、仿真和實(shí)驗(yàn)結(jié)果均體現(xiàn)出了明顯的優(yōu)越性，表明該算法對(duì)非線性生物系統(tǒng)未知參數(shù)的估計(jì)具有較佳的性能。