生物系統參數估計的改進粒子群優化算法

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

(1.河南財政金融學院 計算機與信息技術學院 計算機軟件與理論重點學科，河南 鄭州 450000；2.河南工業大學 糧食信息處理與控制教育部重點實驗室，河南 鄭州 450001)

0 引 言

生物模擬系統的參數估計是一個NP hard問題，具有多模態和病態的特征，被估計的參數通常有多個真實解[1]。梯度優化算法，由于對初始值敏感度高，通常不能得到真實解[2]。啟發式算法，已被證明能夠解決NP hard問題，能夠克服多峰性問題，并得到全局最優解[3]。

生物模擬系統的參數估計問題，有無跡卡爾曼濾波(unscented Kalman filter，UKF)[4,5]、模擬退火(simulated annealing，SA)算法[6]、螢火蟲算法[7]、遺傳算法[8]等。

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

為了驗證算法的有效性，將改進的PSO算法應用于一個仿真的生物模型的參數估計(該生物模型的數據來源于一個E-Coli系統的真實數據)。比較改進PSO算法與迭代無跡卡爾曼濾波器(iterated unscented Kalman filter，IUKF)、SA和PSO算法的RMSE值，結果表明，改進的PSO算法具有更小的RMSE值，從而能夠更精確地估計生物系統仿真模型的參數。

1 數學建模

假定生物系統的非線性動力學模型如下

(1)

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

用來描述生物現象的建模框架有多種可供使用。諸如Lotka-Volterra[9]、S-systems[10]等建模框架中，選擇了S-system來表達式(1)所描述的非線性模型

(2)

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

2 模型求解方法

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

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

UKF是卡爾曼濾波器用于非線性建模的改進框架[11]。在此濾波器中，過程和檢測方程可能是非線性的，因而需要非線性卡爾曼濾波器。UKF算法采用無跡變換作為非線性變換來計算隨機變量的統計量。

UKF算法如下：

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

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

(3)

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

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

(4)

式中：x[k] 和u[k] 為輸入，y[k] 為輸出，w為待估計的具有維數q的未知參數。為了估計參數，定義如下狀態空間

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

(5)

其中，前者為過程方程，r為過程噪聲。后者則是由輸入和測量噪聲e[k] 驅動的測量方程。

UKF算法的收斂速度較慢，濾波器可能不收斂于真實參數的主要原因，是不當的初始值分配，導致了卡爾曼增益的過早飽和。因此，較小的濾波器增益導致參數估計較小的變化，會使算法的收斂變得非常緩慢。

為解決這個問題，開發了IUKF算法。在本算法中，RMSE計算公式如下

(6)

2.2 模擬退火(SA)算法

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

在高溫時的第一次迭代中，與代價函數的非適當值相反，內循環中出現壞解的概率較高。這個特性可防止算法在局部響應中收斂。隨著溫度的降低，在外環的最后迭代中，壞解的概率降低了，算法以較高的概率選擇了最合適的解。

2.3 粒子群算法優化(PSO)

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

(7)

式中：W為慣性權重，通常在0.9～1.2范圍內，c1和c2是設為2的常數，r1和r2為區間[0,1]內的隨機參數。

第i個粒子的新位置更新公式如下

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

(8)

PSO算法是一種快速簡單的方法，適用于非凸NP hard問題，如生物路徑建模。雖然它在搜索解空間上有很好的表現，但在最優解的鄰域中搜索最優解上存在不足。針對這一不足，通常采用一些啟發式算法或將其與其它方法結合來解決。在粒子群算法中嵌入一種新的啟發式算法，提高了算法的尋優能力。在該算法中，在一個序列平臺中使用規范化的PSO開發了一種分解算法，該算法稱為基于分解的PSO算法(particle swarm optimization algorithm base on decomposition，PSOBD)，由兩部分組成。

第一部分，是一個標準的PSO，算法在最優解的鄰域內找到一個解，即全局最優解。為了提高第一部分盡可能接近最優解的能力，在PSO算法的每次迭代中，線性減小式(7)中的慣性權重，最終值為0.1。當粒子在最佳解附近時，使用這種技術，可以在最后的迭代中，提供更精確的解。

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

begin

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

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

-使用PSO算法求第j維的最優值

-更新全局最佳位置

end

end

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

當其它決策變量作為固定的常數值時，PSOBD算法能夠在最佳解的附近，搜索到決策變量的最優值。這種方法對所有決策變量是可重復的。在外部循環中，重復整個過程，以確保避開陷入局部最優。該方法可以產生更高精度的最優解。

3 仿真與實驗

3.1 仿真結果與分析

采用Matlab2016a進行仿真。

通過兩個綜合仿真來說明改進的粒子群算法的性能。在第一個仿真中，假設沒有附加噪聲。在第二個仿真中，考慮了信噪比為20的加性高斯白噪聲。數據通過以下模型生成

(9)

參數如下

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

(10)

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

表2～表5給出了IUKF、SA、PSO和提出的PSOBD算法的仿真結果。在這些表中，給出了每種算法被估參數的1000次仿真平均值和相應的平均RMSE值。PSOBD算法在沒有附加噪聲時RMSE值為0.197，分別比IUKF算法、SA算法和PSO算法下降60.12%、71.16%和49.87%，而在有附加噪聲(SNR=20)時，RMSE值為0.304，分別比IUKF算法、SA算法和PSO算法下降46.95%、64.02%和38.83%。可以看出，PSOBD算法對兩種仿真場景下的未知參數估計，都有較好的效果，該方法的RMSE遠小于其它3種算法。在兩個仿真實驗中，RMSE值降幅達39%～71%。仿真結果表明，與普通粒子群算法相比，PSOBD算法的參數估計能力更強，與IUKF算法和SA算法相比，PSOBD算法具有更精準的參數估計能力。圖2給出了表2～表5中每種算法在有噪聲和無噪聲情況下1000次仿真的RMSE值變化情況。

表1 模擬S系統中選擇的參數

圖1 4種合成狀態的時間分布

相對于SA、IUKF和PSO方法，PSOBD算法在無噪聲和有噪聲場景下的平均百分比分別提升到60.06%和48.72%。

3.2 實驗結果與分析

為了比較算法在實驗場景中的性能，利用了Cad系統中的真實數據集。Cad系統是大腸桿菌條件應力響應模塊之一，僅在低pH和富含賴氨酸的環境下被激活。下面的S系統可以用來模擬這種現象

表2 IUKF算法兩次仿真的估計參數

表3 SA算法兩次仿真的估計參數

表4 PSO算法兩次仿真的估計參數

表5 PSOBD算法兩次仿真的估計參數

圖2 在無噪聲和有噪聲兩種情況下1000次仿真運行的RMSE值

(11)

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

由估計參數得到生成的時間分布和真實時間分布如圖3所示。

由圖3可見，將分解技術應用于一般的粒子群算法中，

表6 4種算法在實際實驗中的RMSE

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

綜合仿真和真實數據集上的實驗，測試結果表明，本算法的均方根誤差比其它對比方法分別平均減少了55.16%和19.62%。

圖3 由估計參數得到的時間分布和真實時間分布

4 討 論

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

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

在實驗中，利用計算機輔助設計系統對一個真實的生物代謝系統進行了建模。該模型與真實場景相似，我們對模型參數應用以上的4種算法進行了估計。實驗結果表明，對于IUKF、SA、PSO和PSOBD算法在CPU為i5 4590，內存8 G條件下，每個算法的運行時間分別為343 s、305 s、287 s和416 s，PSOBD算法運行時間較長，表明PSOBD方法有進一步改進的空間。

5 結束語

提出了一種基于改進粒子群算法的非線性生物模型的未知參數估計算法PSOBD。通過采用分解技術，把PSO算法從多維空間集中到單維空間上，當把其它決策變量作為常數時，PSOBD算法能夠搜索到決策變量的最優值。該算法對所有決策變量是可重復的，在外部循環中，可以避開陷入局部最優，產生更高精度的最優解。將PSOBD算法得到的結果與IUKF算法、SA算法和一般的粒子群優化算法進行了比較，結果表明，與其它3種方法相比，PSOBD算法在速度、仿真和實驗結果均體現出了明顯的優越性，表明該算法對非線性生物系統未知參數的估計具有較佳的性能。