一類生物系統的參數辨識方法及其應用

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

1.渤海大學 數理學院，遼寧 錦州121013

2.渤海大學 數字出版大數據挖掘與治理及呈現技術標準實驗室，遼寧 錦州121013

1 引言

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

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

2 生物系統的數學描述

考慮如下S-型生物系統：

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

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

3 參數辨識問題

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

為了得到生物系統（1）的速率常數和動力階，本文以極小化誤差函數為優化目標，建立了如下參數辨識模型：

參數辨識問題（9）具有如下性質：

性質1參數辨識問題（9）存在最優解p?∈U。

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

4 參數辨識問題的求解方法

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

其中，ηk=tk-tk-1。則參數辨識問題（9）可以轉化為如下非線性規劃問題：

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

非線性規劃問題（11）的求解也可以在Matlab中實現，應用Matlab優化工具箱[11]的fmincon命令即可求得最優解。

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

步驟1構建生物系統（1）的參數辨識優化模型（9）。

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

步驟3應用修正配置公式（10）將生物系統的微分方程（1）近似表示為代數方程，從而將動態優化問題（9）化為非線性規劃問題（11）。

步驟4通過步驟4.1至步驟4.4辨識出生物系統（1）的模型參數。

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

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

步驟5輸出模型參數的辨識值p?。結束參數辨識過程。

5 仿真研究

例1考慮如下生物系統：

生物系統（12）、（13）中各參數的真實值為：α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）可得如下參數辨識優化模型：

表1給出了無噪聲條件下本文方法的參數辨識結果，相應的最優目標函數值為7.446 2×10-8。從表1可以看出，本文方法得到的辨識結果接近參數的真實值。

表1 無噪聲條件下例1的辨識結果

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

表2 例1中誤差e的比較

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

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

圖1 無噪聲條件下例1的仿真結果

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

圖2 噪聲條件下例1的仿真結果

例2考慮如下生物系統：

生物系統（16）～（19）中各參數的真實值為：α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）可得如下參數辨識優化模型：

表3給出了無噪聲條件下本文方法的參數辨識結果，相應的最優目標函數值為3.795 466×10-8。從表3可以看出，本文方法的辨識結果接近參數的真實值。

表3 無噪聲條件下例2的辨識結果

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

表4 例2中誤差e的比較

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

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

圖3 無噪聲條件下例2的仿真結果

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

圖4 噪聲條件下例2的仿真結果

6 應用于甘油生物發酵系統

基于本文方法，研究甘油生物發酵生產1，3-丙二醇系統[12-16]的參數辨識問題。

甘油生物發酵生產1，3-丙二醇的S-型生物系統可表示為：

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

圖5 生物量變化曲線

圖6 甘油濃度變化曲線

圖7 主要產物1，3-丙二醇濃度變化曲線

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

7 結束語

針對S-型生物系統的參數辨識問題，本文給出了以極小化誤差函數為優化目標的參數辨識優化模型，并為其設計了有效的求解方法。與已有方法相比，本文方法能夠獲得更為準確的參數辨識結果。另外，將本文方法應用于甘油生物發酵系統的參數辨識問題中，取得了較好的應用效果。