一類生物系統(tǒng)的參數(shù)辨識方法及其應(yīng)用

徐恭賢，邱 添，賈府生，魏順行

1.渤海大學(xué) 數(shù)理學(xué)院，遼寧 錦州121013

2.渤海大學(xué) 數(shù)字出版大數(shù)據(jù)挖掘與治理及呈現(xiàn)技術(shù)標(biāo)準(zhǔn)實驗室，遼寧 錦州121013

1 引言

隨著生物工程技術(shù)的迅速發(fā)展，越來越多的學(xué)者對S-型生物系統(tǒng)的參數(shù)辨識問題進(jìn)行了深入研究[1-9]。例如，Chen等[1]基于混合變量多目標(biāo)進(jìn)化算法推斷了S-型生物系統(tǒng)的模型參數(shù)；Sarode等[2]利用遺傳算法對S-型生物系統(tǒng)進(jìn)行了反問題研究；Kimura等[3]應(yīng)用分解法把一個原問題分解成若干個子問題，進(jìn)而對每一個子問題進(jìn)行了參數(shù)辨識分析；Tsai等[4]將微分方程轉(zhuǎn)化為代數(shù)方程后應(yīng)用混合差分算法求解參數(shù)辨識的優(yōu)化模型；于超[5]通過五點中心差分算法推斷了S-型生物系統(tǒng)的模型參數(shù)；劉豐[6]、徐恭賢等[7]應(yīng)用三次樣條插值方法研究了S-型生物系統(tǒng)的參數(shù)辨識問題。

本文針對S-型生物系統(tǒng)的參數(shù)辨識問題，基于修正配置和B樣條插值，提出了一種有效的參數(shù)辨識方法。應(yīng)用研究表明，本文方法能夠獲得較好的參數(shù)辨識結(jié)果。

2 生物系統(tǒng)的數(shù)學(xué)描述

考慮如下S-型生物系統(tǒng)：

其中，x=(x1,x2,…,xn)T∈Rn+為生物系統(tǒng)的代謝物濃度向量，t∈[0,T],f=(f1,f2,…,fn)T∈Rn。

式（2）中，αi、βi為速率常數(shù)；gij、hij為動力階；向量p=(αT,βT,GT,HT)T,pL≤p≤pU,p中的α、β、G、H可分別表示為：

3 參數(shù)辨識問題

設(shè)xei(tk)表示第i個代謝物濃度在時刻t=tk時的實驗值；xeimax表示第i個代謝物濃度實驗值的最大值；Ns表示時刻點的個數(shù)；xi(tk)表示第i個代謝物濃度在時刻t=tk時的模型計算值；x?ei(tk)表示第i個代謝物濃度在時刻t=tk時實驗值的速率；x?i(tk)表示第i個代謝物濃度在時刻t=tk時模型計算值的速率；x?eimax表示第i個代謝物濃度速率的最大測量值。

為了得到生物系統(tǒng)（1）的速率常數(shù)和動力階，本文以極小化誤差函數(shù)為優(yōu)化目標(biāo)，建立了如下參數(shù)辨識模型：

參數(shù)辨識問題（9）具有如下性質(zhì)：

性質(zhì)1參數(shù)辨識問題（9）存在最優(yōu)解p?∈U。

證明 設(shè)U={p|pL≤p≤pU,p∈，顯然U?是緊集。又因為x0∈Rn+,p∈U且fi∈C1(Rn+×U)，故由常微分方程的解存在唯一性定理可知，對?p∈U，生物系統(tǒng)（1）存在唯一解，且該解關(guān)于p∈U是連續(xù)的，所以J關(guān)于x是連續(xù)泛函，故參數(shù)辨識問題（9）存在最優(yōu)解p?∈U。

4 參數(shù)辨識問題的求解方法

為了求得參數(shù)辨識問題（9）的最優(yōu)解，利用修正配置法將微分方程近似表示為如下代數(shù)方程：

其中，ηk=tk-tk-1。則參數(shù)辨識問題（9）可以轉(zhuǎn)化為如下非線性規(guī)劃問題：

為了求解非線性規(guī)劃問題（11），本文應(yīng)用B樣條插值方法估計實驗值的速率B樣條插值是一類常見的樣條函數(shù)[10]，一般選擇階次為4或5的B樣條能夠得到較好的插值效果。關(guān)于B樣條插值的基本原理及其求解過程可以參見文獻(xiàn)[10]。B樣條插值方法的應(yīng)用可以在Matlab中完成，采用Matlab樣條插值工具箱[11]的spapi命令即可計算實驗值的速率

非線性規(guī)劃問題（11）的求解也可以在Matlab中實現(xiàn)，應(yīng)用Matlab優(yōu)化工具箱[11]的fmincon命令即可求得最優(yōu)解。

綜合前文所述，本文的參數(shù)辨識方法可描述如下：

步驟1構(gòu)建生物系統(tǒng)（1）的參數(shù)辨識優(yōu)化模型（9）。

步驟2應(yīng)用B樣條插值（比如Matlab樣條插值工具箱的spapi命令）估計實驗值的速率

步驟3應(yīng)用修正配置公式（10）將生物系統(tǒng)的微分方程（1）近似表示為代數(shù)方程，從而將動態(tài)優(yōu)化問題（9）化為非線性規(guī)劃問題（11）。

步驟4通過步驟4.1至步驟4.4辨識出生物系統(tǒng)（1）的模型參數(shù)。

步驟4.1取初始參數(shù)值p(0)∈U，ε＞0。令迭代次數(shù)r=0。

步驟4.2在算法的第r(r≥1)次迭代，應(yīng)用非線性優(yōu)化方法的序列二次規(guī)劃算法求解非線性規(guī)劃問題（11）。設(shè)動力階參數(shù)gij和hij的最優(yōu)值分別為和。

步驟5輸出模型參數(shù)的辨識值p?。結(jié)束參數(shù)辨識過程。

5 仿真研究

例1考慮如下生物系統(tǒng)：

生物系統(tǒng)（12）、（13）中各參數(shù)的真實值為：α1=3；α2=1；β1=1；β2=1；g11=0；g12=-2；g21=0.5；g22=1；h11=0.5；h12=1；h21=0；h22=0.5。本例中，取x(0)=(0.2,0.5)T，T=8，ηk=0.1。

由式（9）可得如下參數(shù)辨識優(yōu)化模型：

表1給出了無噪聲條件下本文方法的參數(shù)辨識結(jié)果，相應(yīng)的最優(yōu)目標(biāo)函數(shù)值為7.446 2×10-8。從表1可以看出，本文方法得到的辨識結(jié)果接近參數(shù)的真實值。

表1 無噪聲條件下例1的辨識結(jié)果

表2給出了本文方法與已有方法[5-6]的結(jié)果比較。誤差函數(shù)e表示參數(shù)辨識值與真實值之差的絕對值之和，即：

表2 例1中誤差e的比較

通過表2可以看出，本文方法獲得了較低的誤差e值，分別比文獻(xiàn)[5-6]方法降低了32.86%和0.21%，說明本文方法更為精確。

將表1的參數(shù)辨識結(jié)果代入生物系統(tǒng)（12）、（13）中，得到如圖1所示的代謝物濃度隨時間變化曲線。從圖1可以看出，各代謝物濃度的計算值與實驗值基本一致。

圖1 無噪聲條件下例1的仿真結(jié)果

為考察本文方法在有噪聲條件下的性能，在代謝物濃度的實驗值xei(tk)中加入5%的高斯白噪聲，參數(shù)辨識后生物系統(tǒng)的仿真結(jié)果如圖2所示。從圖2中可以看出，在噪聲情況下，本文方法仍能得到較好的參數(shù)辨識結(jié)果。

圖2 噪聲條件下例1的仿真結(jié)果

例2考慮如下生物系統(tǒng)：

生物系統(tǒng)（16）～（19）中各參數(shù)的真實值為：α1=12；α2=8；α3=3；α4=2；β1=10；β2=3；β3=5；β4=6；g11=0；g12=0；g13=-0.8；g14=0；g21=0.5；g22=0；g23=0；g24=0；g31=0；g32=0.75；g33=0；g34=0；g41=0.5；g42=0；g43=0；g44=0；h11=0.5；h12=0；h13=0；h14=0；h21=0；h22=0.75；h23=0；h24=0；h31=0；h32=0；h33=0.5；h34=0.2；h41=0；h42=0；h43=0；h44=0.8。本例中，取x(0)=(10,1,2,3)T，T=5，ηk=0.1。

由式（9）可得如下參數(shù)辨識優(yōu)化模型：

表3給出了無噪聲條件下本文方法的參數(shù)辨識結(jié)果，相應(yīng)的最優(yōu)目標(biāo)函數(shù)值為3.795 466×10-8。從表3可以看出，本文方法的辨識結(jié)果接近參數(shù)的真實值。

表3 無噪聲條件下例2的辨識結(jié)果

表4給出了本文方法與已有方法[5-6]的結(jié)果比較。誤差函數(shù)e表示參數(shù)辨識值與真實值之差的絕對值之和，即：

表4 例2中誤差e的比較

通過表4可以看出，本文方法獲得了與文獻(xiàn)[6]方法相同的誤差e值，但比文獻(xiàn)[5]方法降低了1.93%。

將表3的參數(shù)辨識結(jié)果代入生物系統(tǒng)（16）～（19）中，得到如圖3所示的代謝物濃度隨時間變化曲線。從圖3可以看出，各代謝物濃度的計算值與實驗值基本一致。

圖3 無噪聲條件下例2的仿真結(jié)果

為考察本文方法在有噪聲條件下的性能，在代謝物濃度的實驗值xei(tk)中加入5%的高斯白噪聲，參數(shù)辨識后生物系統(tǒng)的仿真結(jié)果如圖4所示。從圖4中可以看出，在噪聲情況下，本文方法仍能得到較好的參數(shù)辨識結(jié)果。

圖4 噪聲條件下例2的仿真結(jié)果

6 應(yīng)用于甘油生物發(fā)酵系統(tǒng)

基于本文方法，研究甘油生物發(fā)酵生產(chǎn)1，3-丙二醇系統(tǒng)[12-16]的參數(shù)辨識問題。

甘油生物發(fā)酵生產(chǎn)1，3-丙二醇的S-型生物系統(tǒng)可表示為：

其中，x1表示生物量的濃度(g?L-1)，x2表示甘油的濃度(mmol?L-1)，x3表示1，3-丙二醇的濃度(mmol?L-1)，x4表示乙酸的濃度(mmol?L-1)，x5表示乙醇的濃度(mmol?L-1)。式（22）～（26）的參數(shù)辨識優(yōu)化問題為：

圖5 生物量變化曲線

圖6 甘油濃度變化曲線

圖7 主要產(chǎn)物1，3-丙二醇濃度變化曲線

本例中，取x0=(0.1905,400.043,0,0,0)T,T=6。圖5～圖7分別表示甘油生物發(fā)酵過程中生物量、甘油、主要產(chǎn)物1，3-丙二醇濃度隨時間的變化曲線，其中離散點表示實驗的觀測數(shù)據(jù)。由擬合曲線可以看出，應(yīng)用本文方法得到的S-系統(tǒng)能較好地描述甘油生物發(fā)酵系統(tǒng)。

7 結(jié)束語

針對S-型生物系統(tǒng)的參數(shù)辨識問題，本文給出了以極小化誤差函數(shù)為優(yōu)化目標(biāo)的參數(shù)辨識優(yōu)化模型，并為其設(shè)計了有效的求解方法。與已有方法相比，本文方法能夠獲得更為準(zhǔn)確的參數(shù)辨識結(jié)果。另外，將本文方法應(yīng)用于甘油生物發(fā)酵系統(tǒng)的參數(shù)辨識問題中，取得了較好的應(yīng)用效果。