利用符號(hào)計(jì)算方法研究生物系統(tǒng)全時(shí)滯穩(wěn)定性

柳君，牛薇

(北京航空航天大學(xué) 中法工程師學(xué)院，北京100191)

時(shí)滯微分系統(tǒng)的穩(wěn)定性研究在理論和應(yīng)用上都有其重要的意義，特別地，從控制理論的角度看，生物系統(tǒng)全時(shí)滯穩(wěn)定即表明該系統(tǒng)對(duì)于時(shí)滯具有很好的魯棒性和可靠性.

長(zhǎng)期以來(lái)，人們一直致力于尋找時(shí)滯微分系統(tǒng)全時(shí)滯穩(wěn)定的代數(shù)判據(jù)，已取得了不少進(jìn)展.秦元?jiǎng)自谖墨I(xiàn)［1］中第1 次將單滯后多維系統(tǒng)的全時(shí)滯穩(wěn)定判據(jù)由超越形式的檢驗(yàn)轉(zhuǎn)化為代數(shù)形式的檢驗(yàn)，并對(duì)n =1，2(n 為非線性時(shí)滯微分系統(tǒng)的維數(shù))的情形具體給出了判定法則，雖然此方法對(duì)高維情形處理起來(lái)并不容易.在文獻(xiàn)［2］中，秦元?jiǎng)缀陀嵩楦鶕?jù)Newton 遞推公式，對(duì)單時(shí)滯n=2 的情況進(jìn)行了深入的研究.文獻(xiàn)［3］討論了一類線性定常時(shí)滯系統(tǒng)全時(shí)滯漸進(jìn)穩(wěn)定的充分代數(shù)判據(jù).文獻(xiàn)［4-5］研究了單時(shí)滯線性系統(tǒng)漸進(jìn)穩(wěn)定的代數(shù)判據(jù);在文獻(xiàn)［5］的基礎(chǔ)上，文獻(xiàn)［6］討論了多時(shí)滯線性系統(tǒng)漸進(jìn)穩(wěn)定性的代數(shù)判據(jù).文獻(xiàn)［1-6］中的方法均針對(duì)線性系統(tǒng)，并均使用傳統(tǒng)的數(shù)學(xué)推導(dǎo)方法.模擬生物問(wèn)題的動(dòng)力系統(tǒng)通常是非線性的，較為復(fù)雜，利用傳統(tǒng)的數(shù)學(xué)方法很難分析其穩(wěn)定性，需要借助更先進(jìn)的計(jì)算方法和工具.隨著計(jì)算機(jī)科學(xué)與技術(shù)的發(fā)展，以精確計(jì)算為特點(diǎn)的符號(hào)計(jì)算逐漸成熟和完善，效率也逐步提高，成為數(shù)值計(jì)算的一種強(qiáng)有力的替代.

針對(duì)含參數(shù)的非線性生物系統(tǒng)，本文給出時(shí)滯微分系統(tǒng)全時(shí)滯穩(wěn)定性的代數(shù)判據(jù)，研究如何利用符號(hào)計(jì)算方法分析多維生物系統(tǒng)正平衡點(diǎn)的漸進(jìn)穩(wěn)定性.文獻(xiàn)［7-8］提出了利用代數(shù)方法分析生物系統(tǒng)穩(wěn)定性和分岔等問(wèn)題的算法化方法，并且給出了軟件實(shí)現(xiàn)及實(shí)驗(yàn)結(jié)果.不同于傳統(tǒng)的數(shù)學(xué)推導(dǎo)，符號(hào)計(jì)算使得求解生物系統(tǒng)的穩(wěn)定性和分岔更加程序化和自動(dòng)化.本文在作者工作的基礎(chǔ)上，進(jìn)一步研究如何利用符號(hào)計(jì)算方法，如文獻(xiàn)［9］中的三角分解、文獻(xiàn)［10-11］中的Gr?bner 基和文獻(xiàn)［12］中的實(shí)解分類等方法，程序化地分析生物系統(tǒng)全時(shí)滯穩(wěn)定性.

1 時(shí)滯微分系統(tǒng)正平衡點(diǎn)的穩(wěn)定性

考慮如下n 維非線性多時(shí)滯微分系統(tǒng):

式中:P1，P2，…，Pn，Q1≠0，Q2≠0，…，Qn≠0 為多項(xiàng)式;u = (u1，u2，…，um)為參數(shù);x = (x1(t)，x2(t)，…，xn(t))∈Rn為變?cè)?τ=(τ，2τ，…，kτ)為系統(tǒng)的時(shí)滯，τ∈R+=［0，∞)為時(shí)滯;n、m、k為正整數(shù).

為了計(jì)算式(1)的正平衡點(diǎn)，令τ=0，可得

得出的解即為平衡點(diǎn)x*.在求出系統(tǒng)(1)的正平衡點(diǎn)之后，需要討論其平衡點(diǎn)x*處的穩(wěn)定性.由于該系統(tǒng)是非線性形式，首先需進(jìn)行線性化.據(jù)文獻(xiàn)［13］，線性化之后的系統(tǒng)零解的漸進(jìn)穩(wěn)定性與非線性系統(tǒng)正平衡點(diǎn)的漸進(jìn)穩(wěn)定性一致.

令y =x －x*，代入式(1)中，x*是系統(tǒng)的正平衡點(diǎn)，滿足式(2)，可得線性化之后的系統(tǒng):

式中:A，Ak∈Rn×n為矩陣.

式(3)的特征方程為

式中:I 為n×n 階單位矩陣;λ 為特征值.若特征方程的根對(duì)?τ∈R+均具有負(fù)實(shí)部，那么其零解漸進(jìn)穩(wěn)定，稱式(3)系統(tǒng)全時(shí)滯穩(wěn)定.因此，式(1)系統(tǒng)正平衡點(diǎn)的漸進(jìn)穩(wěn)定性等價(jià)于式(3)系統(tǒng)零解的全時(shí)滯穩(wěn)性.

2 時(shí)滯系統(tǒng)全時(shí)滯穩(wěn)定的代數(shù)判據(jù)

時(shí)滯微分式(3)系統(tǒng)的全時(shí)滯穩(wěn)定性，即多項(xiàng)式Δ(λ，τ)對(duì)?τ∈R+均Hurwitz 穩(wěn)定.據(jù)文獻(xiàn)［14］，有

引理1 式(3)系統(tǒng)全時(shí)滯穩(wěn)定的充分必要條件為

2)對(duì)?τ ＞0 及任意實(shí)數(shù)y，都有

成立.其中，x*為式(2)求出的平衡點(diǎn);i 為虛數(shù)單位.

對(duì)于條件1):令多項(xiàng)式M 為

在此基礎(chǔ)上，定義M 的Hurwitz 矩陣:

其中:當(dāng)i ＞n 時(shí)，ai=0.H 的順序主子式Δ1，Δ2，…，Δn為M 的Hurwitz 行列式.根據(jù)文獻(xiàn)［15］Routh-Hurwitz 判據(jù)，多項(xiàng)式M 是穩(wěn)定的，當(dāng)且僅當(dāng)下列條件成立:

Δ1＞0，Δ2＞0，…，Δn－1＞0，an＞0

對(duì)于條件2):對(duì)?τ ＞0 及任意y，有

成立.據(jù)文獻(xiàn)［14］，令－yτ =θ，則其他等價(jià)于:對(duì)?θ∈［0，2π］＞0 及任意實(shí)數(shù)y∈R－{0}，有

成立.作分式線性變換

條件2)等價(jià)于:對(duì)?z∈R 及任意實(shí)數(shù)y∈R－{0}，都有

成立.在此，記

由h(u，z，y)=0 分離實(shí)部和虛部可得

式中:f 和g 為關(guān)于z、y 的實(shí)系數(shù)含參數(shù)二元多項(xiàng)式，z∈R，y∈R－{0};Re 和Im 為實(shí)部和虛部.條件2)等價(jià)于式(4)無(wú)實(shí)根.由此，可得非線性多時(shí)滯微分系統(tǒng)全時(shí)滯穩(wěn)定的充要條件.

引理2 式(1)系統(tǒng)正平衡點(diǎn)全時(shí)滯穩(wěn)定的充分必要條件是式(5)系統(tǒng)有正實(shí)根且式(6)系統(tǒng)無(wú)實(shí)根:

3 代數(shù)方法分析

第2 節(jié)描述了如何將時(shí)滯系統(tǒng)的全時(shí)滯穩(wěn)定性問(wèn)題化為代數(shù)問(wèn)題，這里將研究如何將代數(shù)方法應(yīng)用到求解這些代數(shù)問(wèn)題上.

步驟1 構(gòu)造半代數(shù)系統(tǒng).

假設(shè)實(shí)際問(wèn)題由式(1)系統(tǒng)來(lái)模擬，通過(guò)計(jì)算可得到含多項(xiàng)式方程及不等式的式(5)系統(tǒng)和式(6)系統(tǒng).變?cè)獂 及參數(shù)u 可能需滿足一些實(shí)際問(wèn)題中的附加約束:

式中:s，t≥n.把約束條件式(7)加到式(5)系統(tǒng)和式(6)系統(tǒng)中，可得兩個(gè)含等式及不等式的系統(tǒng)，稱之為半代數(shù)系統(tǒng).一般地，設(shè)半代數(shù)系統(tǒng)Ψ形如:

步驟2 計(jì)算三角列.

通過(guò)三角列或Gr?bner 基方法，可以把多項(xiàng)式集E={E1，E2，…，Es}三角化，得到三角列Tk.若參數(shù)u 出現(xiàn)，轉(zhuǎn)至步驟4.

步驟3 參數(shù)不出現(xiàn).

如果參數(shù)u 不出現(xiàn)，可以用長(zhǎng)度可無(wú)限小的有理區(qū)間來(lái)隔離每個(gè)Tk的實(shí)零點(diǎn).令F = {N1，N2，…，Nt，H1，H2，…，Hn}為不等式多項(xiàng)式集.然后F 中的多項(xiàng)式在每個(gè)實(shí)零點(diǎn)上的符號(hào)可以通過(guò)計(jì)算它們?cè)谟欣韰^(qū)間的端點(diǎn)的值來(lái)確定.由此得到用有理區(qū)間隔離的半代數(shù)系統(tǒng)Ψ 的實(shí)解.

步驟4 實(shí)解分類.

假設(shè)參數(shù)u 出現(xiàn)，令F 對(duì)于每個(gè)三角列Tk，可利用F 來(lái)計(jì)算一個(gè)關(guān)于u 的代數(shù)簇V，使得其將實(shí)參數(shù)空間Rm分為有限多個(gè)胞腔，滿足在每個(gè)胞腔上，Tk的實(shí)零點(diǎn)個(gè)數(shù)以及F 在這些零點(diǎn)處的符號(hào)不變.從每個(gè)胞腔中取一個(gè)有理樣本點(diǎn)，代入Tk消去參數(shù)，由步驟3 可計(jì)算Tk的實(shí)零點(diǎn)以及F 在該樣本點(diǎn)的符號(hào).

步驟5 參數(shù)條件最后，根據(jù)多項(xiàng)式在Ψ 具有指定個(gè)數(shù)的實(shí)解的胞腔中樣本點(diǎn)處的符號(hào)，建立參數(shù)u 所需要滿足的條件.

在文獻(xiàn)［16］實(shí)解分類軟件DISCOVERER 的基礎(chǔ)上，算法步驟(見(jiàn)圖1)已在MAPLE 中實(shí)現(xiàn).

圖1 代數(shù)方法分析算法步驟Fig.1 Algorithmic steps using algebraic approaches

4 應(yīng)用實(shí)例

4.1 實(shí)例驗(yàn)證

文獻(xiàn)［13］中介紹了一種重要的生物模型:Lotka-Volterra 捕食-食餌系統(tǒng).為了驗(yàn)證本文中代數(shù)方法的可行性，使用這個(gè)簡(jiǎn)單的捕食-食餌系統(tǒng)演示利用符號(hào)計(jì)算方法分析全時(shí)滯穩(wěn)定性的流程.在此，考慮以下單時(shí)滯的捕食-食餌系統(tǒng):

式中:x(t)＞0、y(t)＞0 為食餌、捕食者的種群密度;r1＞0、r2＞0 為食餌及捕食者的內(nèi)稟增長(zhǎng)率;aii(i=1，2)＞0 為兩種群密度作用的種群內(nèi)作用系數(shù);aij(i≠j)＞0 為兩種群相互作用的種群間作用系數(shù);τ≥0 為捕食者的追捕時(shí)間.接下來(lái)計(jì)算該系統(tǒng)在平衡點(diǎn)處穩(wěn)定的充要條件.

第1 步 分析正平衡點(diǎn)及線性化.

令τ=0，那么平衡點(diǎn)問(wèn)題可化為

記由此求得的正平衡點(diǎn)x*=(x*，y*).

考慮2 ×2 階Jacobian 矩陣

式(8)系統(tǒng)可寫為矩陣形式:

這里G 是非線性項(xiàng)，線性部分為

式中:

第2 步 分析條件1).

第2 節(jié)引理1 中的條件1)是det［λI －A －A1］=0 均具有負(fù)實(shí)根，在此，

M=det［λI － A － A1］=a0λ2+ a1λ + a2

式中:

由多項(xiàng)式系數(shù)a0、a1、a2組成的Hurwitz 判據(jù)

故而，可得半代數(shù)系統(tǒng)

利用MAPLE 中實(shí)現(xiàn)的第3 節(jié)中的算法求解該半代數(shù)系統(tǒng)可得

R1=r1a21－ r2a11＞0

第3 步 分析條件2).

無(wú)實(shí)根.

分離h(y，z)的實(shí)部和虛部可得半代數(shù)系統(tǒng)

無(wú)實(shí)根.

求解該半代數(shù)系統(tǒng)可得

R2=a21a12－ a11a22≤0

滿足時(shí)，系統(tǒng)無(wú)實(shí)根.

綜上可得，當(dāng)參數(shù)滿足

時(shí)式(8)系統(tǒng)在平衡點(diǎn)處全時(shí)滯穩(wěn)定.

4.2 無(wú)選擇性的捕食-食餌系統(tǒng)

文獻(xiàn)［17］中的無(wú)選擇收獲的Lotka-Volterra捕食-食餌系統(tǒng)也是一種重要的數(shù)學(xué)生物模型.考慮以下二維系統(tǒng):

式中:r(1 －x/k)(r，k ＞0)為無(wú)捕食者時(shí)食餌的增長(zhǎng)率;βx/(1 + ax)(β，a ＞0)為捕食者的響應(yīng)函數(shù);d ＞0為捕食者的死亡率;α ＞0 為轉(zhuǎn)換系數(shù);qx和Ey(q，E ＞0)分別為食餌和捕食者的收獲率.與第4.1 節(jié)的方法類似可以解決式(9)系統(tǒng)的全時(shí)滯穩(wěn)定的問(wèn)題.

第1 步 利用線性化方法和Hurwitz 判據(jù)可得半代數(shù)系統(tǒng)，即條件1):

求解該半代數(shù)系統(tǒng)可得

第2 步 分析條件2).

無(wú)實(shí)根.

分離h(y，z)的實(shí)部和虛部可得半代數(shù)系統(tǒng)

無(wú)實(shí)根.

f(y，z)和g(y，z)都是含有參數(shù)和變量的多項(xiàng)式.表1 中是f(y，z)和g(y，z)的項(xiàng)數(shù)和最高次數(shù)，該問(wèn)題使用傳統(tǒng)的數(shù)學(xué)計(jì)算方法很難得出結(jié)果，而利用符號(hào)計(jì)算方法在幾秒之內(nèi)，可得

R4= － ark + aqk － r≥0

滿足時(shí)，系統(tǒng)無(wú)實(shí)根.

表1 f(y，z)和g(y，z)的項(xiàng)數(shù)和最高次數(shù)Table 1 Number of terms and higher degree of f(y，z)and g(y，z)

所以，當(dāng)參數(shù)滿足

時(shí)，式(9)系統(tǒng)在平衡點(diǎn)處全時(shí)滯穩(wěn)定.

4.3 SIR 傳染病模型

SIR 傳染病模型是一個(gè)重要的生物模型，Cooke 在文獻(xiàn)［18］中提出了時(shí)滯SIR 傳染病模型，并指出t 時(shí)刻的傳染能力為βS(t)I(t －τ')，其中，β ＞0 為每天每個(gè)感染者接觸的人數(shù)，τ'≥0為病毒在被感染者體內(nèi)的作用時(shí)間.文獻(xiàn)［19］中時(shí)滯SIR 傳染病系統(tǒng)為

式中:μ ＞0 為死亡率;r ＞0 為日恢復(fù)速率.

第1 步 構(gòu)造及分析條件1).

計(jì)算可得

R1=β － μ － r ＞0

第2 步 分析條件2).

分離h(y，z)的實(shí)部和虛部可得半代數(shù)系統(tǒng)

其中:f(y，z)有30 項(xiàng)，最高次數(shù)是7;g(y，z)有29項(xiàng)，最高次數(shù)是7.通過(guò)計(jì)算可得參數(shù)取任意值時(shí)該系統(tǒng)均無(wú)實(shí)根.

因此，當(dāng)參數(shù)滿足R1= β － μ － r ＞0 時(shí)，式(10)系統(tǒng)在平衡點(diǎn)處全時(shí)滯穩(wěn)定.

5 結(jié) 論

本文在Gr?bner 基、三角化分解和實(shí)解分類等符號(hào)計(jì)算原理基礎(chǔ)上提出了一種新的驗(yàn)證生物系統(tǒng)全時(shí)滯穩(wěn)定性的算法.

1)經(jīng)實(shí)驗(yàn)驗(yàn)證表明該算法可實(shí)現(xiàn)較為優(yōu)異的計(jì)算性能，例如計(jì)算含有62 項(xiàng)的多項(xiàng)式所用時(shí)間僅僅為幾秒，這是傳統(tǒng)的數(shù)學(xué)計(jì)算所達(dá)不到的.

2)此外，仍在進(jìn)行任意多時(shí)滯微分系統(tǒng)的全時(shí)滯穩(wěn)定性分析的算法研究及實(shí)驗(yàn).

References)

［1］ 秦元?jiǎng)?有時(shí)滯的系統(tǒng)的無(wú)條件穩(wěn)定性［J］.數(shù)學(xué)學(xué)報(bào)，1960，10(1):125-142.

Chin Y S.Unconditional stability of systems with time lags［J］.Acta Mathematica Sinica，1960，10(1):125-142(in Chinese).

［2］ 秦元?jiǎng)?，俞元?一類時(shí)滯微分系統(tǒng)無(wú)條件穩(wěn)定的條件［J］.控制理論與應(yīng)用，1984，1(1):23-35.

Chin Y S，Yu Y H.Unconditional stability conditions for a class of differential systems with time delay［J］.Journal of Control Theory and Applications，1984，1(1):23-35(in Chinese).

［3］ 周超順，鄧聚龍.線性定常時(shí)滯系統(tǒng)全時(shí)滯漸近穩(wěn)定的充分代數(shù)判據(jù)［J］.自動(dòng)化學(xué)報(bào)，1990，16(1):62-65.

Zhou C S，Deng J L.A sufficient algebra criteria for stability of linear constant time-delay system［J］.Acta Automatica Sinica，1990，16(1):62-65(in Chinese).

［4］ Cao D Q，He P，Ge Y M.Simple algebraic criteria for stability of neutral delay-differential systems［J］.Journal of the Franklin Institute，2005，342(3):311-320.

［5］ Hu G D，Hu G D，Cahlon B.Algebraic criteria for stability of linear neutral systems with a single delay［J］.Journal of Computational and Applied Mathematics，2001，135(1):125-133.

［6］ He P，Cao D Q.Algebraic stability criteria of linear neutral systems with multiple time delays［J］.Applied Mathematics and Computation，2004，155(3):643-653.

［7］ Niu W，Wang D.Algebraic analysis of bifurcation and limit cycles for biological systems［M］.Berlin，Heidelberg:Springer，2008:156-171.

［8］ Niu W，Wang D.Algebraic approaches to stability analysis of biological systems［J］.Mathematics in Computer Science，2008，1(3):507-539.

［9］ Wang D.Elimination methods［M］.Berlin，Heidelberg:Springer，2001:193-224.

［10］ Buchberger B.Gr?bner bases:An algorithmic method in polynomial ideal theory［J］.Multidimensional Systems Theory，1985:184-232.

［11］ Faugère J C.A new efficient algorithm for computing Gr?bner bases (F4)［J］.Journal of Pure and Applied Algebra，1999，139(1):61-88.

［12］ Yang L，Xia B.Real solution classification for parametric semialgebraic systems［C］∥Algorithmic Algebra and Logic［S.l.:s.n.］，2005:281-289.

［13］ 宋永利，韓茂安，魏俊杰.多時(shí)滯捕食-食餌系統(tǒng)正平衡點(diǎn)的穩(wěn)定性及全局Hopf 分支［J］.數(shù)學(xué)年刊:A 輯，2005，25(6):783-790.

Song Y L，Han M A，Wei J J.Stability and global Hopf bifurcation for a predator-prey model with two delays［J］.Chinese Annals of Mathematics:Series A，2005，25(6):783-790(in Chinese).

［14］ Bhattacharyya S P，Chapellat H，Keel L H.Robust control-the parametric approach［M］.New York:Prentice Hall PTR，1995:446-472.

［15］ Lancaster P，Tismenetsky M.The theory of matrices:With applications［M］.Pittsburgh:Academic Press，1985:89-103.

［16］ Xia B.DISCOVERER:A tool for solving semi-algebraic systems［J］.ACM Communications in Computer Algebra，2007，41(3):102-103.

［17］ Kar T K，Pahari U K.Non-selective harvesting in prey-predator models with delay［J］.Communications in Nonlinear Science and Numerical Simulation，2006，11(4):499-509.

［18］ Cooke K L.Stability analysis for a vector disease model［J］.Journal of Mathmatics，1979，9(1):31-42.

［19］ Meng X，Chen L，Wu B.A delay SIR epidemic model with pulse vaccination and incubation times［J］.Nonlinear Analysis:Real World Applications，2010，11(1):88-98.